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SUMMARY 


This  thesis  presents  a  robust  method  for  estimating  the  relaxations  of  a  metallic  object 
from  its  electromagnetic  induction  (EMI)  response.  The  EMI  response  of  a  metallic  object 
can  be  accurately  modeled  by  a  sum  of  real  decaying  exponentials.  However,  it  is  difficult  to 
obtain  the  model  parameters  from  measurements  when  the  number  of  exponentials  in  the 
sum  is  unknown  or  the  terms  are  strongly  correlated.  Traditionally,  the  relaxation  constants 
are  estimated  by  nonlinear  iterative  search  that  often  leads  to  unsatisfactory  results. 

An  effective  EMI  modeling  technique  is  developed  by  first  linearizing  the  problem 
through  enumeration  and  then  solving  the  linearized  model  using  a  sparsity-regularized 
minimization.  This  approach  overcomes  several  long-standing  challenges  in  EMI  signal 
modeling,  including  finding  the  unknown  model  order  as  well  as  handling  the  ill-posed 
nature  of  the  problem.  The  resulting  algorithm  does  not  require  a  good  initial  guess  to 
converge  to  a  satisfactory  solution. 

This  new  modeling  technique  is  extended  to  incorporate  multiple  measurements  in  a 
single  parameter  estimation  step.  More  accurate  estimates  are  obtained  by  exploiting  an 
invariance  property  of  the  EMI  response,  which  states  that  the  relaxation  frequencies  do  not 
change  for  different  locations  and  orientations  of  a  metallic  object.  Using  tests  on  synthetic 
data  and  laboratory  measurement  of  known  targets,  the  proposed  multiple-measurement 
method  is  shown  to  provide  accurate  and  stable  estimates  of  the  model  parameters. 

The  ability  to  estimate  the  relaxation  constants  of  targets  enables  more  robust  subsur¬ 
face  target  discrimination  using  the  relaxations.  A  simple  relaxation-based  subsurface  target 
detection  algorithm  is  developed  to  demonstrate  the  potential  of  the  estimated  relaxations. 
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CHAPTER  I 


INTRODUCTION 


The  landmine  crisis  remains  today  as  landmines  continue  to  maim  or  kill  civilians  everyday 
worldwide.  The  International  Campaign  to  Ban  Landmines  reported  that  in  the  year  of 
2009,  landmines  and  explosive  remnants  of  war  caused  about  4000  casualties  worldwide,  of 
which  over  60%  are  civilians  [1]  and  more  than  30%  are  children.  Much  effort  and  research 
has  been  invested  in  remediating  landmines  with  one  of  the  primary  tasks  being  the  detection 
of  the  landmine  itself.  However,  landmine  detection  can  suffer  from  a  high  false-alarm  rate 
as  the  detectors  also  detect  other  metallic  non-mine  objects  like  gun  shells,  metal  cans,  and 
shrapnel.  Therefore,  it  is  of  strong  interest  to  discriminate  between  landmines  and  metallic 
non-mine  objects. 

Recent  research  has  shown  that  broadband  electromagnetic  induction  (EMI)  sensors 
are  capable  of  discriminating  between  certain  types  of  targets  [2-7].  Target  discrimination 
using  broadband  EMI  sensors  is  possible  because  the  EMI  response  of  a  target  is  strongly 
related  to  the  target’s  physical  size,  shape,  orientation,  and  composition. 

EMI  sensors  work  by  illuminating  a  target  with  a  time- varying  magnetic  field,  and  then 
detecting  the  scattered  magnetic  field  which  is  generated  by  the  eddy  currents  induced 
on  the  target.  The  broadband  EMI  sensors  measure  the  scattered  field  at  a  finite  set  of 
frequencies  or  over  a  short  interval  measurement  times.  The  measured  response  can  be  fitted 
to  a  parametric  model,  and  discrimination  of  targets  is  performed  based  on  the  fitted  model 
parameters.  This  research  is  performed  in  the  context  of  frequency-domain  broadband  EMI 
sensors,  which  have  the  advantage  of  avoiding  undesired  signal  sources,  such  as  power  line 
harmonics. 

The  main  contribution  of  this  work  is  a  new  effective  EMI  modeling  technique  developed 
by  first  linearizing  through  enumeration  and  then  solving  the  linearized  model  using  a 
sparsity-regularized  minimization.  Furthermore,  this  approach  is  extended  to  incorporate 
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multiple  measurements  simultaneously,  utilizing  physical  properties  of  the  EMI  signal  in 
the  modeling  process,  from  which  more  accurate  estimates  are  obtained.  This  approach 
overcomes  several  long  standing  challenges  in  EMI  signal  modeling,  including  the  unknown 
model  order  as  well  as  the  ill-posed  condition  of  the  problem. 

1.1  EMI  Measurement  System 

This  research  is  conducted  as  part  of  the  development  of  a  practical  field  measurement 
system  to  find  buried  targets.  Shown  in  Fig.  1  is  a  measurement  cart  equipped  with  three 
frequency-domain  broadband  EMI  sensors  [8].  In  the  held,  this  cart  is  mobile,  and  when  it 
is  driven  near  a  target,  the  EMI  frequency  response  of  the  target  is  acquired  sequentially,  as 
illustrated  in  Fig.  2.  As  the  cart  moves,  multiple  EMI  responses  of  a  target  are  measured  at 
different  relative  sensor-to-target  positions  and  locations.  In  addition  to  the  target  response, 
the  frequency  response  of  the  soil  is  also  measured  because  soil  is  always  present  in  the  held 
and  interacts  with  the  EMI  system.  Therefore,  the  total  response  can  be  described  as 

Total  response  =  Target  response  +  Soil  response  +  System  noise.  (1) 

The  system  noise  is  primarily  due  to  the  thermal  noise  introduced  in  the  preampliher  of  the 
system,  which  can  be  modeled  as  a  Gaussian  white  noise  [9]. 


Figure  1:  A  broadband  EMI  system  mounted  on  a  cart. 
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Figure  2:  Sequential  data  acquisition  of  the  EMI  system. 


To  detect  and  discriminate  subsurface  targets  in  the  scenario  described  above,  a  detec¬ 
tion  framework,  shown  in  Fig.  3,  is  developed  in  this  research.  Given  a  measured  response, 
the  detection  framework  first  passes  the  response  through  a  soil  prescreener  to  determine 
if  only  soil  is  present.  If  a  soil-only  decision  cannot  be  made,  then  a  metallic  object  might 
be  present,  and  this  response  is  then  passed  to  an  estimator  where  the  response  is  fitted  to 
an  EMI  model  and  the  spectrum  of  the  object  is  estimated.  Lastly,  a  classifier  using  the 
estimated  spectrum  as  a  feature  vector  decides  whether  this  response  is  from  a  landmine  or 
clutter. 

* 

Figure  3:  Block  diagram  for  the  a  system  that  classifies  target  from  measured  EMI  re¬ 
sponses. 
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This  research  focuses  on  developing  techniques  to  model  the  EMI  response,  or,  equiva¬ 
lently,  to  estimate  its  underlying  spectrum.  The  background  for  EMI  modeling  is  presented 
in  Section  1.2.  To  demonstrate  target  discrimination  using  the  modeled  response,  classi¬ 
fiers  and  a  soil  prescreener  are  also  developed.  The  background  for  target  discrimination  is 
presented  in  Section  1.3. 

1 . 2  EMI  Models 

Several  EMI  models  have  been  proposed,  including  ones  for  canonical  targets  [10, 11]  and 
ones  for  targets  of  more  general  shape  [12-14],  In  general,  these  models  are  instances  of  a 
Fredholm  integral  equation  of  the  first  kind: 

H(u)  =  g0  +  gA  f°°  dr,  (2) 

Jo  1  +  Jwr 
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where  H (u/)  is  the  frequency  response,  r  the  relaxation  time,  and  go  and  are  constants; 
G(t)  is  the  distribution  function,  or  the  spectrum. 

Depending  on  characteristics  of  the  distribution  function,  the  EMI  models  assume  a 
distribution  of  relaxations  that  is  either  continuous  or  discrete.  These  two  types  of  spectra 
are  examined  in  the  following  sections,  where  it  is  also  shown  that  the  discrete  spectrum 
is  appealing  for  the  work  in  EMI  discrimination  because  of  its  physical  basis  as  well  as  its 
invariant  properties. 

1.2.1  Continuous  Distribution  of  Relaxations 


Miller  et  al.  proposed  a  four-parameter  model  [14],  which  is  also  identified  with  the  Cole- 
Cole  model  used  in  polymer  science  [15]: 


Hcciu)  =  go  + 


a  a 


(3) 


1  +  (jUT0)Q  ’ 

where  To  and  a  are  model  parameters.  It  is  known  from  polymer  science  that  this  model 
has  a  continuous  distribution  of  relaxation  times  (DRT) 

1  sin(a7r) 


Gcc(r)  = 


27 r  cosh(a  log(r/ro))  +  cos(a7r) 


(4) 


The  Cole-Cole  model  assumes  a  DRT  that  is  symmetric  in  logr.  However,  not  all 
targets  have  a  symmetric  DRT,  e.g.  a  sphere  [16].  More  complex  DRTs  can  be  described  by 
more  complex  models,  such  as  the  Cole-Davidson  or  the  Havriliak-Negami  model  [17,18]. 
Examples  of  these  DRT  models  can  be  found  in  Fig.  6.  While  the  Cole-Cole  model  may 
not  exactly  describe  a  target,  it  has  been  shown  to  be  practical  for  target  discrimination  in 
the  context  of  EMI  [6,7]. 


1.2.2  Discrete  Spectrum  of  Relaxations 


Several  researchers  have  provided  a  theoretical  basis  for  modeling  the  EMI  response  of  a 
metallic  object  with  an  underlying  discrete  distribution  of  relaxations  [12,13].  In  this  case, 
the  unit-step  time  response  h(t)  is  a  discrete  sum  of  damped  real  exponentials: 


h(t)  =  dou(t)  —  ^^dke  t^Tku(t), 
k= 1 


(5) 
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where  Tfe  >  0  are  the  real  relaxation  time  constants,  d &  the  real-valued  coefficients,  and  u(t) 
the  unit  step  function. 


In  the  frequency  domain,  the  EMI  response  can  be  written  as 

K 


h(u)  =  c0 + 


1  -  jWCfc 
^  ^  1  +  jVCfc  ’ 


(6) 


where  ( k  =  7fc-1  are  the  relaxation  frequencies,  Co  the  shift,  Ck  the  real  spectral  amplitudes, 
and  K  the  model  order.  Comparing  (5)  to  (6),  the  model  parameters  are  related  by  do  = 
J2k= ock  and  dk  =  2 Cfc,  where  do  =  H( 0)  is  due  to  the  dc  magnetization  of  the  target.  In 
this  work,  the  parameter  set  S  =  {(Cfc>  c&)  :  k  =  1 . . .  K\  is  called  the  Discrete  Spectrum  of 
Relaxation  Frequencies  (DSRF)  or  simply  the  spectrum;  each  pair  (£*,,  c& )  is  defined  as  one 
relaxation. 

To  visualize  the  time-domain  and  the  frequency-domain  model  and  their  relationship, 
consider  a  response  with  S  =  { ( 104,  0.5),  (105, 0.5)}  and  Co  =  —1.  The  time-domain  response 
is  shown  in  Fig.  4,  where  h(t)  is  a  sum  of  two  decaying  exponentials.  The  frequency  response 
H(lu)  measured  at  21  frequencies  is  shown  in  Fig.  5,  where  the  real  and  imaginary  parts 
of  H(lu)  are  shown  separately  as  well  as  together  on  an  Argand  diagram  (complex  plane), 
parameterized  by  ui.  The  plots  of  H(u)  in  Fig.  5  illustrate  an  important  fact  about  modeling 
the  EMI  response  in  which  only  a  limited  range  of  the  response  is  observable.  This  constraint 
influences  how  accurately  the  true  model  of  an  EMI  response  can  be  recovered.  It  is  difficult 
to  recover  relaxations  that  are  far  outside  of  the  observable  range  because  their  contributions 
to  the  response  may  not  be  observed  or  are  indistinguishable  from  the  constant  cq. 


Figure  4:  An  example  unit-step  time  response  h(t). 
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Figure  5:  An  example  frequency  response  H(oj)  sampled  at  21  frequencies  over  300-90 kHz.  The  dotted  lines  trace  out  the  full  range 
of  the  response.  The  solid  lines  trace  out  the  measured  response,  i.e.,  the  observable  part  for  a  band-limited  EMI  system.  The  solid 
circles  indicate  the  samples  measured  at  21  distinct  frequencies,  (a)  The  real  part  and  (b)  the  imaginary  part  of  H(oj).  (c)  The  complex 
response,  parameterized  by  w,  plotted  on  an  Argand  diagram  (complex  plane),  (d)  The  DSRF,  where  each  relaxation  is  represented  by 
a  stem:  the  stem  location  is  the  relaxation  frequency  Cfc  and  the  stem  height  is  the  spectral  amplitude  Ck-  The  dash-dotted  lines  in  (d) 
indicate  the  observable  range. 


A  detailed  derivation  of  the  relation  between  (5)  and  (6),  as  well  as  other  forms  of  the 
EMI  model,  can  be  found  in  Appendix  A.  The  significance  of  choosing  the  particular  form 
in  (6)  is  that  when  this  expression  is  linearized,  as  in  the  proposed  method,  the  linearized 
model  results  in  a  matrix  with  uniform-norm  columns.  All  columns  having  the  same  norm 
translates  to  an  unbiased  estimation  of  relaxation  frequencies,  as  explained  in  Appendix  A. 

The  measured  EMI  response  is  proportional  to  the  product  of  the  transmitted  and 
received  magnetic  fields  and  the  magnetic  polarizability  tensor  of  the  target  being  measured 
by  the  EMI  sensor  (Appendix  A).  The  magnetic  polarizability  tensor  of  several  canonical 
targets  can  be  calculated  analytically,  and  these  formulas  show  how  the  model  parameters 
are  related  to  the  target’s  physical  properties  such  as  conductivity,  permeability,  shape, 
size,  and  orientation  [12,16,19,20]. 

One  result  of  this  analysis  is  that  the  relaxation  frequencies  are  invariant  with  re¬ 
spect  to  the  target’s  position  and  orientation  relative  to  the  EMI  sensor;  only  the  spectral 
amplitudes  change  with  orientation  and  position.  Therefore,  the  DSRF,  i.e. ,  the  ([^’s  and 
Cfc’s,  are  valuable  for  target  discrimination,  and  even  allow  estimation  of  the  target  depth 
and  orientation. 

1.2.3  Existing  Estimation  Methods  and  Challenges 

It  is  difficult  to  find  a  meaningful  low-order  model  for  an  EMI  frequency  response.  This  is 
because  the  modeling  process  is  an  ill-posed  problem  -  a  small  change  in  H{oS)  can  lead  to 
a  large  change  in  the  DSRF  and  different  DSRFs  can  generate  very  similar  or  even  identical 
responses.  To  illustrate  the  ill-posed  nature  of  the  modeling  problem,  shown  in  Fig.  6  are 
curve  fits  for  the  same  target  response  using  different  spectral  models.  It  is  clear  that  all 
four  models  (spectra)  generate  a  response  very  close  to  the  measurement.  However,  the 
four  spectra  are  actually  very  different,  which  illustrates  the  ill-posed  nature  of  the  inverse 
problem. 

Several  spectrum  modeling  techniques  have  been  developed  in  the  past  few  decades  for 
both  continuous  and  discrete  spectra.  For  continuous  spectra,  modeling  techniques  such 
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as  Tikhonov  regularization  [21],  a  nonparametric  Bayesian  approach  [22],  and  a  Monte- 
Carlo  method  [23]  have  been  proposed.  In  particular,  for  the  Cole-Cole  model,  existing 
methods  include  nonlinear  iterative  search,  multifold  least-squares  estimation  [24],  and  a 
gradient-lookup-table  approach  [7]. 
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Figure  6:  A  target  response  (left)  fitted  to  four  different  spectral  models  (right):  Cole- 
Cole,  Havriliak-Negami,  a  general  DRT,  and  the  DSRF,  plus  a  standard  least-squares  fit. 
All  models  fit  well  but  have  very  different  distribution  function  G(t).  The  DRT  fit  is 
obtained  using  the  Tikhonov  regularization. 


In  the  discrete  case,  the  modeling  problem  is  equivalent  to  estimating  the  parameters 
of  a  sum  of  exponentials,  which  is  a  classic  problem  that  arises  in  many  fields.  Several 


notable  estimation  techniques  include  iterative  nonlinear  least-squares  fitting,  the  matrix 
pencil  method  [25,26],  modified  Prony’s  methods  [27,28]  ,  and  complex-curve  fitting  [29]. 
While  the  listed  methods  can  be  robust  in  other  fields  of  science,  in  the  case  of  EMI  spectral 
modeling,  a  major  challenge  is  that  the  poles  are  required  to  be  real-valued ,  but  most  of 
these  model-fitting  methods  can  return  complex  poles.  In  addition,  often  these  methods  do 
not  perform  well  when  three  or  more  relaxations  are  present.  The  goodness  of  estimation 
strongly  depends  on  a  good  guess  of  the  model  order,  and  is  also  very  sensitive  to  the  initial 
guess  for  the  model  parameters. 

Riggs  et  al.  proposed  a  differential  correction  technique  that  returns  only  real  relaxations, 
but  the  convergence  still  strongly  depends  on  a  good  initial  guess  of  the  model  order.  When 
the  guess  of  K  is  poor,  the  techniques  fails  to  converge  [5].  Estimating  the  unknown  model 
order  is  the  primary  source  of  DSRF  estimation  difficulties. 

In  addition,  as  a  consequence  of  the  ill-posed  problem,  it  is  possible  to  converge  well  to 
a  small  residual,  but  obtain  a  fit  that  is,  in  fact,  far  from  the  correct  solution.  This  is  a 
common  problem  in  the  DSRF  estimation  to  obtain  “good  fits  but  poor  estimates,”  where 
the  fitting  residual  can  be  minimized  but  the  estimated  DSRF  is  far  off.  The  least-squares 
fit  in  Fig.  6  is  a  good  example  of  this  scenario,  where  the  fitting  residual  is  the  smallest  of  all 
methods,  but  the  estimated  spectrum  is  noise-like  with  oscillatory,  self-canceling  relaxations. 

The  nonlinear  relationship  between  H(uj )  and  (  often  requires  a  nonlinear  iterative 
search  to  fit  a  response  to  the  model.  Some  challenges  of  the  search  are  that  1)  the  iteration 
can  settle  in  local  minima  that  are  far  from  the  correct  solution,  2)  the  solution  is  sensitive 
to  the  initial  guess,  3)  the  returned  solution  may  be  complex  and  lack  physical  meaning.  It 
is  possible  that  an  estimate  is  only  a  good  numerical  fit,  but  has  no  physical  significance  [30]. 

Das  and  McFee  provided  a  detailed  examination  of  these  challenges  [30].  In  summary, 
the  challenges  of  DSRF  estimation  are  1)  estimating  the  model  order  K ,  2)  not  converging, 
3)  good  numerical  fit  but  poor  estimate,  and  4)  returning  complex  estimates. 

For  these  reasons,  a  simplified  DSRF  model  with  only  one  or  two  relaxations  (K  =  1  or 
2)  is  often  assumed  [2,31,32].  The  goal  of  this  work  is  to  accurately  estimate  parameters 
from  the  full  model  without  constraining  or  assuming  K. 
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1.2.4  Nonnegative  Least-Squares  (NNLSQ) 


The  initial  effort  in  this  research  suggests  a  constrained  linear  method  that  robustly  esti¬ 
mates  the  DSRFs  without  requiring  prior  knowledge  of  the  model  order  and  bypasses  the 
nonlinear  relationship  between  £  and  H(u>)  [33,34],  This  is  achieved  through  enumeration 
to  obtain  a  linear  problem,  and  then  minimization  of  the  squared  error  with  a  nonnega¬ 
tive  constraint  on  the  spectral  amplitudes.  The  nonnegative  constraint  effectively  reduces 
the  size  of  the  feasible  solution  set  and  eliminates  oscillatory  and  nonphysical  solutions. 
For  example,  the  oscillatory  least-squares  solution  in  Fig.  6  is  rejected  by  the  nonnegative 
constraint. 

While  the  NNLSQ  algorithm  performs  well,  the  method  presumes  nonnegative  spectra 
for  the  targets,  i.e.,  Ck  >  0.  This  assumption  is,  in  fact,  valid  for  most  targets  measured 
by  our  system  [8,33],  but  it  is  possible  that  some  targets  can  have  a  spectrum  with  both 
positive  and  negative  relaxations.  A  goal  of  this  work  is  to  develop  a  more  general  method 
that  remove  the  nonnegative  constraint. 

1.2.5  Proposed  Estimation  Methods 

To  robustly  estimate  the  DSRF  without  the  nonnegative  constraint  but  still  require  no 
prior  knowledge  of  K  and  return  only  real- valued  solutions,  as  in  the  NNLSQ,  a  sparsity  - 
promoting  DSRF  estimation  method  is  formulated.  The  sparsity  approach  follows  the  same 
sampling  and  linearization  setup  as  in  the  NNLSQ,  but  replaces  the  nonnegative  constraint 
with  a  sparsity-regularization  term.  The  sparsity  regularization  is  shown  to  deliver  robust 
estimates  without  prior  knowledge  of  K  and  returns  only  real-valued  estimates. 

This  work  further  extends  the  DSRF  estimation  from  the  single-measurement  case  to 
estimation  using  multiple  measurements.  In  the  laboratory  and  the  field,  multiple  mea¬ 
surements  are  often  available  for  a  given  target  of  interest.  Because  different  measurements 
from  the  same  target  share  the  same  relaxation  frequencies,  this  invariance  property  can  be 
exploited  to  increase  the  estimation  accuracy  when  multiple  measurements  are  available. 

The  multiple  measurements  are  utilized  simultaneously  by  recasting  the  EMI  estima¬ 
tion  problem  into  the  problem  of  multiple-measurement  vectors  (MMV)  or  jointly-sparse 
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vectors,  where  vectors  sharing  the  same  underlying  structure  are  recovered.  The  model 
parameters  can  then  be  estimated  through  iterative  reweighting  algorithms.  The  proposed 
MMV  method  is  tested  against  synthetic,  laboratory,  and  field  data,  and  is  demonstrated 
to  deliver  accurate  and  stable  estimates. 

To  account  for  the  soil  response  in  held  measurements,  the  proposed  DSRF  estimation 
method  is  also  extended  to  model  the  soil.  A  soil  model  is  proposed  based  on  the  observation 
of  a  large  number  of  soil  responses,  of  which  a  subset  is  shown  in  Fig.  7.  From  the  soil 
measurements,  it  is  observed  that  the  frequency  responses  follow  a  simple  log-linear  trend. 
That  is, 

#soiiH  =  7  (\n—  +  (7) 

\  coo  2  J 

where  7  and  wo  are  model  parameters.  This  model  has  a  linear  real  part  with  respect  to 
the  log-frequency  and  has  a  constant  imaginary  part,  which  agrees  with  the  observed  soil 
measurements.  Using  this  soil  model,  an  augmented  model  is  developed  to  consider  both 
the  DSRF  and  the  soil  simultaneously.  It  is  shown  that  the  augmented  model  can  also  be 
solved  by  the  same  proposed  sparsity-promoting  DSRF  estimation  method. 


x  10"7 


Figure  7:  The  frequency  response  of  soil  samples  plotted  on  an  Argand  diagram.  The 
real  and  imaginary  part  are  shown  in  Fig.  21.  The  soil  responses  are  collected  at  various 
locations  in  a  testing  field. 
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1.3  EMI-based,  Target  Discrimination 

1.3.1  Existing  Discrimination  Methods 

Several  EMI-based  landmine  discrimination  techniques  have  been  developed  in  the  past 
decade.  Both  Bayesian  and  non-Bayesian  classification  as  well  as  the  continuous  and  discrete 
spectral  model  have  been  employed.  Fails  et  al.  [6]  and  Ramachandran  et  al.  [7]  both 
demonstrated  success  in  detecting  mines  using  the  k- nearest-neighbor  (kNN)  algorithm 
based  on  the  EMI  model  developed  by  Miller  et  al.,  i.e.,  the  Cole-Cole  model  [14,15].  Using 
the  Cole-Cole  model,  the  target  feature  is  a  vector  in  M4.  In  [7],  the  modeling  error  is  also 
included  in  the  feature  vector.  The  kNN  labels  an  unknown  target  based  on  the  closest 
targets  in  a  training  set. 

Collins  et  al.  [2]  and  Gao  et  al.  [4]  considered  the  likelihood  ratio  test  applied  to  both 
time  and  frequency  domain  data.  In  [2],  it  is  considered  using  only  the  most  prominent 
relaxation,  i.e.,  K  =  1  in  (5).  In  [4],  a  body-of-revolution  model  is  considered. 

Because  the  Cole-Cole  spectrum  and  the  most  prominent  relaxation  are  dependent  on 
the  orientation  and  location  of  a  target,  the  above  discrimination  approaches  can  suffer 
from  orientation  variations  of  buried  targets.  The  use  of  the  DSRF  takes  into  account  this 
variation  and  is  considered  by  Riggs  et  al.  [5]  and  is  also  considered  in  this  work. 

Riggs  et  al.  suggested  identifying  targets  using  a  library  of  relaxation  frequencies  ob¬ 
tained  from  training  targets  [5].  This  approach  takes  into  account  the  orientation  and 
position  variation.  The  method  uses  a  generalized  likelihood  ratio  test,  which  results  in  a 
simple  minimum  least-squares-error  decision  rule.  In  a  sense,  this  test  is  a  kNN  approach 
that  considers  only  the  closest  neighbor. 

1.3.2  Proposed  Discrimination  Method 

In  this  work,  a  discrimination  based  on  the  estimated  DSRF  is  considered.  Figure  8  demon¬ 
strates  the  consistency  of  the  estimated  DSRF  across  instances  of  two  different  types  of 
targets.  For  each  target  type,  the  EMI  responses  are  collect  from  independent  physical  in¬ 
stances  buried  at  different  locations  and  depths.  It  is  seen  that  instances  from  the  same  type 
share  a  consistent  DSRF,  while  the  two  types  clearly  have  distinct  DSRFs.  This  observation 
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agrees  with  Baum’s  theory  [12]  and  is  the  motivation  of  DSRF-based  target  discrimination. 


(b) 
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Figure  8:  Estimated  DSRF  of  landmines.  The  spectral  amplitude  Ck  is  represented  by  the 
intensity:  the  darker  the  gray,  the  larger  the  amplitude,  (a)  Seven  low-metal  content,  non¬ 
magnetic,  moderate  EMI  response  antipersonnel  mines,  (b)  Eight  medium-metal  content, 
magnetic,  strong  EMI  response  antipersonnel  mines. 


Two  classifiers  are  examined  for  the  proposed  DSRF  discrimination:  kNN  and  the  sup¬ 
port  vector  machine  (SVM).  The  kNN  serves  as  a  simple,  proof-of-concept  DSRF-based 
classification.  The  SVM  is  a  more  practical  approach  as  the  computation  cost  is  lower  for 
real-time  applications. 

To  build  a  complete  system  for  discriminating  targets,  a  model-based  soil  prescreener  is 
also  developed  to  detect  the  presence  of  metallic  objects.  The  prescreener  is  based  on  the 
soil  model  (7),  and  the  prescreener  turns  out  to  be  a  simple  least-squares  minimizer,  which 
can  be  computed  efficiently.  Combining  the  prescreener  and  the  DSRF-based  classifier,  the 
proposed  discrimination  technique  is  shown  to  be  competitive  to  existing  methods. 
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1.4  Outline 


The  outline  of  this  thesis  is  as  follows:  the  DSRF  estimation  method  based  on  a  single 
measurement  vector  (SMV)  utilizing  sparsity  is  presented  in  Chapter  II;  this  SMV  method 
is  extended  to  multiple  measurements  in  Chapter  III;  the  application  of  DSRF  in  target 
discrimination  is  demonstrated  in  Chapter  IV,  where  the  soil  prescreener  is  also  introduced; 
lastly,  a  conclusion  is  provided  in  Chapter  V. 
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CHAPTER  II 


SPECTRUM  ESTIMATION  FROM  SINGLE  MEASUREMENTS 


To  robustly  estimate  the  DSRF  without  the  nonnegative  constraint  from  a  given  EMI  fre¬ 
quency  response,  a  sparsity-promoting  estimation  method  is  developed.  The  estimation 
problem  is  first  recast  into  a  linear  dictionary  selection  problem  by  densely  sampling  the 
relaxation  frequency  space.  The  linearized  problem  is  then  solved  by  minimizing  the  least- 
squares  error  regularized  with  sparsity.  Nonconvex  iterative  sparsity-promoting  optimiza¬ 
tion  is  employed  to  achieve  robust  performance.  An  interpolation  is  performed  after  the 
optimization  step  to  address  the  issue  if  an  actual  relaxation  frequency  is  split  between  two 
sampled  relaxation  frequencies  [33]. 

The  proposed  sparsity-promoting  DSRF  estimation  method  is  shown  to  be  robust  through 
tests  on  synthetic,  laboratory,  and  field  measurements.  The  method  gives  accurate  estimates 
and  requires  no  a  priori  knowledge  of  the  model  order  and  always  returns  real  model  param¬ 
eters.  An  extension  to  the  Earth  Mover’s  Distance  (EMD)  [35]  is  devised  to  quantify  the 
estimation  accuracy.  A  simulation  framework  is  also  proposed  to  find  the  best  regularization 
parameter  used  in  the  sparsity-promoting  optimization. 

The  soil  response  is  also  considered  when  modeling  the  field  measurements.  In  the 
field,  a  target  response  is  always  measured  in  the  presence  of  soil,  so  the  total  response 
is  composed  of  the  target  response,  the  soil  response  and  the  system  noise.  To  take  into 
account  the  effect  of  the  soil,  a  soil  model  is  proposed  based  on  collected  soil  measurements. 
An  augmented  model  is  formulated  to  consider  both  the  DSRF  and  soil  concurrently.  This 
augmented  model  is  also  solved  by  the  sparsity-regularized  least-squares. 

The  DSRF  of  many  field  targets  from  different  types  are  examined  using  the  proposed 
estimation  method,  and  it  is  observed  that  targets  of  the  same  type  have  a  similar  DSRF. 
This  observation  affirms  the  idea  to  use  DSRF  as  a  signature  for  target  discrimination,  and 
this  is  explored  in  Chapter  4. 
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2.1  Method  Formulation 

When  the  target  response  is  measured  at  N  distinct  frequencies,  (6)  can  be  written  in 
matrix  form: 

Hg/ Ci  1-.M/C2  i-jun/c  k 

i+j^i/Ci  1+.M/C2  ' ' ' 

l-i^2/Cl  l-j^2/C2  l-jq'2/Cx 

l+j^/Cl  l+j^2/C2  '  '  '  l+i^2/Cx 

1-j^iv/Cl  1—  j^N/C.2  1—  juN /C,K 

1+i^jv/Ci  l+iwiv/C2  ■  '  '  I+juinKk 

'  V - 

z 

b  =  Zc,  (8) 

where  cum;n  =  ui  <  <  ■  ■  ■  <  ujn  =  wmax,  b  £  Cw  is  the  observation  vector,  c  €  M^+1 

the  spectral  amplitude  vector  augmented  by  the  shift  Co,  and  Z  £  (j^Nx^+i)  a  matrix 
containing  information  about  the  relaxation  frequencies  C-  The  dimension  of  the  matrix  Z 
is  dependent  on  the  number  of  relaxations  present  in  the  spectrum,  i.e. ,  the  model  order. 
In  the  case  of  a  simple  thin  wire  circular  loop,  there  is  only  one  relaxation,  so  Z  has  two 
columns;  the  first  column  is  always  one  to  account  for  cq. 

2.1.1  Sampling  The  Relaxation  Frequency  Domain 

To  linearize  the  DSRF  estimation  problem,  a  large  set  of  sampled  relaxation  frequencies 
{Cm  :  m  =  1 ...  M},  M  3>  K,  are  generated  by  sampling  within  a  range  of  possible  or  ob¬ 
servable  relaxation  frequencies  [Cmim  Cmax]-  There  are  several  ways  to  sample  the  relaxation 
frequency  domain.  The  goal  is  to  sample  densely  enough  to  ensure  an  actual  relaxation 
frequency  Cfc  will  be  close  to  some  sampled  relaxation  frequencies  Cm- 

Because  the  EMI  data  are  measured  at  a  wide  range  of  frequencies,  the  observable 
relaxation  frequencies  also  fall  in  a  wide  range  of  frequency.  To  sample  in  a  high  dynamic 
range,  the  relaxation  frequencies  are  uniformly  logarithmically  spaced,  i.e., 

log  Cm+i  =  log  Cm  +  A,  (9) 

where  A  is  a  fixed  spacing  between  two  C  in  log-space. 


H(u  i)  1 

H{u2)  _  1 

H(ujn)  1 
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This  sampling  scheme  turns  out  to  be  similar  to  the  one  proposed  by  Austin  et  al.  [36], 
where  the  sampling  is  nonuniform  based  on  the  Fisher  information: 

Cm+l  =  Cm  T  OL 

The  parameter  a  controls  the  number  of  samples  given  a  frequency  range.  The  two  sampling 
schemes  produce  similar  sample  points  (Fig.  9).  Uniform  sampling  is  used  in  this  work 
because  it  is  simpler. 


Figure  9:  Samples  generated  from  uniform  log-<C  sampling  and  non-uniform  sampling  based 
on  the  Fisher  information. 


The  expression  (10)  comes  from  [36] 

Cm+l  =  Cm  +  al(b(C m))~1/2, 


(11) 


where  X  is  the  Fisher  information  and  b( (m)  £  is  the  frequency  response  due  to  one 
relaxation  frequency  Cm  measured  at  N  frequencies.  Assuming  the  measurement  noise  is 
Gaussian,  then  [36] 

x(b(C))  ~  MCfMO,  (12) 


where  Jb(C)  is  the  Jacobian  of  6(C)  and  the  superscript  H  denotes  the  Hermitian  transpose. 
The  Jacobian  is 


where 


MC)  = 


Q  0  Q 

^b(C-un)  —„b{ C;w2)  ...  ^ 
d(  d(  d( 


6(C;  ujn) 


(13) 


~b((',  w„) 
d( 


d  1  -jun/C 
d(  l+jun/C 
2jujn 

(c  +  JWn) 


(14) 
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With  some  simple  algebra,  it  can  be  shown  that 

MCmfMU)  =  Y 


n=  1 
N 


2  jujn 


(< 


+  jun 


UJr 


=  V4 

W  + 

Substituting  (15)  into  (11),  (10)  is  obtained. 

2.1.2  Linearizing  The  Problem 

Using  the  sampled  relaxation  frequencies,  (8)  can  be  rewritten  as 


ff(wi) 

1 

l-JWl/Cl 

1-JW1/C2 

1— jwi /Cm 

1+JWl/Cl 

1+JW1/C2 

1+jwi/Cm 

1 

1-JW2/C1 

1-JW2/C2 

1— JW2 /Cm 

= 

1+JW2/C1 

1+1^2/C2 

1+1^2 /Cm 

H(un)_ 

1 

1-jwjv/Cl 

1  j^N  /  C2 

1-/<Wv/Cm 

S  >s 

1+JWjv/Cl 

l+i^jv/C2 

1+JWAt/Cm. 

co 

ci 

C2 

CM 


+e 


(15) 


or 


b  =  Ax  +  e, 


(16) 


where  cm  G  M  are  the  spectral  amplitude  estimates  corresponding  to  each  Cm,  A  G  CArx^A/+1^) 
the  overcomplete  dictionary,  and  e  G  the  modeling  error.  The  solution  vector  x  G  RM+1 
is  a  weighted  selection  vector  containing  the  shift  estimator  co  followed  by  the  spectral  am¬ 
plitude  estimators.  The  matrices  Z  and  A  differ  in  that  Z  is  constructed  with  the  actual 
relaxation  frequencies  while  A  is  the  dictionary  constructed  with  sampled  relaxation  fre¬ 
quencies. 

Because  M  3>  K,  a  good  solution  for  x  will  have  many  zero  elements,  i.e.,  x  will  be 
sparse.  The  intent  is  to  use  x  to  pick  out  the  (m  that  are  close  to  the  actual  relaxation 
frequencies  Cfc-  By  assigning  a  zero  value  to  cm  that  correspond  to  the  Cm  that  are  not 
near  any  Cfc,  the  corresponding  Cm  can  be  made  ineffective.  When  the  zeros  are  properly 
assigned  and  Cfc  properly  assigned  to  cm,  (8)  and  (16)  becomes  almost  identical  and  the 
error  is  small.  In  this  case,  x  is  sparse,  having  only  a  few  nonzero  entries. 
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2.1.3  Sparsity-promoting  Minimization 


While  it  is  true  that  a  good  solution  for  x  produces  a  small  error  in  (16),  directly  minimizing 
the  error  often  yields  solutions  far  from  the  correct  answer.  The  challenge  in  obtaining  the 
correct  x  is  that  M  is  much  greater  than  N,  so  the  system  in  (16)  is  underdetermined  and 
there  is  not  an  unique  x  to  minimize  the  error: 

argmin  ||b  —  Acc|||.  (17) 

X 

Any  vector  in  the  null  space  of  A  can  be  added  to  x  without  changing  the  error.  There 
are  many  ways  to  select  a  least-squares  (LSQ)  solution.  The  Moore-Penrose  pseudoinverse 
picks  the  LSQ  solution  that  has  the  smallest  i 2  norm.  One  can  also  compute  a  LSQ  solution 
with  the  fewest  nonzero  components.  However,  neither  of  these  LSQ  solutions  produces  the 
correct  spectrum,  as  demonstrated  in  Fig.  6.  One  reason  is  that  the  error  in  (16)  will 
not  be  zero  due  to  modeling  error  because  it  is  unlikely  that  all  Q  overlap  with  some  ( m . 
Details  about  existing  techniques  and  the  difficulties  of  solving  such  a  system  can  be  found 
in  [30,37,38]. 

Knowing  that  the  solution  vector  x  should  be  sparse,  the  LSQ  optimization  (17)  can 
be  regularized  to  promote  sparsity.  Sparsity-regularized  least  squares  turn  out  to  be  an 
effective  way  to  obtain  good  solutions  for  x. 

There  are  several  ways  to  perform  sparsity-regularized  least  squares,  including  basis 
pursuit  (7 1  -minimization) ,  matching  pursuit,  iteratively  reweighted  Q -minimization  (IRL1), 
and  iteratively  reweighted  least  squares  (IRLS  or  FOCUSS).  It  has  been  shown  that  the 
latter  two  nonconvex  iterative  reweighted  algorithms,  IRL1  and  FOCUSS,  produce  solutions 
that  are  more  sparse  compared  to  the  non- reweighted  ones  [39].  In  addition,  IRL1  and 
FOCUSS  are  equivalent  under  certain  circumstances.  To  simplify  the  discussion,  only  IR.L1 
is  considered  here.  The  relationship  between  IRL1  and  FOCUSS  is  discussed  in  detail  in 
Chapter  3,  where  multiple  measurements  are  considered. 

To  estimate  the  DSRF  using  sparsity-regularized  least  squares,  the  following  optimiza¬ 
tion  problem  is  solved: 

argmin  \\b'  —  A/ad],  +  A||a:||0  ,0<p<l  (18) 

X  ^ 
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Xt(A)  SHe(6) 

where  A  =  ,  b  = 

3m(A)  3m(b) 

and  A  is  the  regularization  parameter.  Separating  the  real  and  imaginary  parts  in  A  makes 
the  whole  system  real.  Ideally,  in  the  optimal  x,  only  those  cm  with  corresponding  Qm  that 
are  near  a  true  Ck  will  be  nonzero,  and  they  will  take  on  the  correct  spectral  amplitudes 
Cfc.  It  follows  that  a  DSRF  can  then  be  deduced  from  the  nonzero  estimated  cm  and  their 
corresponding  (m. 

The  ip- regularized  least  squares  solution  for  p  <  1  can  be  approximated  by  the  itera¬ 
tively  reweighted  l\  algorithm  proposed  by  Candes  et  al.  [40] .  The  weights  are  updated  as 
suggested  in  [39],  and  the  e-regularization  technique  is  adopted  from  the  same  paper.  In 
summary,  (18)  is  approximated  by  (see  also  [41]): 

Algorithm  1:  Approximated  ^-regularized  least  squares 

Input:  A' ,  b! ,  p,  A,  x° 

1  Xn  4—  X° 

2  for  k  4—  0  to  —8  step  —1  do 

3  e  4-  10fc 

4  repeat 

5  Xn_1 4-  Xn 

6  <^-(|c?"1|  +  e)P-1 

7  xn  arg  min  1 1  b'  —  A'x  \  \  \  +  A  \  Ci  \ 

8  until  \\xn  —  xn~l\\2  <  A/e/100 

9  return  xn 

The  t\  minimization  problem  in  step  7  is  solved  by  ll_ls,  a  Matlab  optimizer  proposed 
by  Kim  et  al.  [42],  It  has  been  observed  that  normalizing  the  input  data  b,  as  well  as  the 
columns  of  A' ,  to  have  an  uniform  £ 2  norm  increases  the  accuracy  of  estimation.  This  is 
why  the  EMI  model  in  the  form  in  (6)  is  chosen  over  other  forms  -  it  produces  a  uniform- 
column-norm  dictionary;  details  are  provided  in  Appendix  A.  While  it  is  often  suggested  to 
initialize  a;0  using  the  least-squares  solution  [39] ,  it  is  observed  that  setting  entries  of  x°  to 
all  ones  also  seems  to  be  effective  and  converges  faster.  The  nonzero  entries  of  x  selected 
by  (18)  along  with  the  corresponding  are  the  relaxations  needed  in  the  estimated  DSRF, 

£  =  {(&, 4)  :  k  =  l...K}. 
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2.2  Estimation  Results  With  Synthetic  Data 


The  proposed  estimation  method  is  tested  against  synthetic  and  laboratory  data  to  show 
its  functionality,  accuracy,  and  stability.  The  hardware  system  used  in  laboratory  and 
field  measurements  is  a  wideband  EMI  sensor  operating  at  21  frequencies  approximately 
logarithmically  distributed  over  the  range  300  Hz-90  kHz,  which  is  a  span  of  2.5  decades  [8]. 
The  synthetic  data  is  generated  in  accordance  with  the  hardware  specification. 

The  range  of  C,  for  estimation  is  chosen  such  that  log(Cmin)  and  log(£max)  are  2.45  and 
6.62,  respectively,  i.e.,  a  span  of  4.17  decades.  All  estimations  are  performed  with  M  =  100, 
and  all  presented  spectra  are  normalized  such  that  ^  \a\  =  1.  Spectral  amplitudes  less  than 
10~5  are  not  displayed.  All  frequency  responses  are  normalized  such  that  | j £»| 1 2  =  1.  Unless 
specified,  p  =  0.5  is  chosen  as  a  representative  case.  In  assessing  the  signal  strength,  the 
signal-to-noise  ratio  (SNR)  is  used.  The  signal  power  is  computed  by  YliLi  |-H(wi)|2/./V. 
The  noise  power  in  synthesized  data  is  equal  to  the  variance  of  the  noise. 

The  regularization  parameter  A  is  chosen  based  on  the  method  described  in  Section  2.5. 
Results  presented  in  this  section  may  achieve  higher  accuracy  with  a  more  sophisticated  A 
selection  rule.  The  purpose  of  this  section  is  to  demonstrate  the  usability  of  the  proposed 
algorithm  using  a  simple  A  selection  rule. 

Notation :  £  and  c  are  the  true/theoretical  relaxation  frequencies  and  spectral  amplitudes; 
(  and  c  are  the  estimates. 

2.2.1  Dissimilarity  Measure  Between  Two  DSRFs 

To  evaluate  the  goodness  of  estimation,  it  is  necessary  to  define  a  measure  to  quantify  the 
similarity  between  two  DSRFs.  It  is  not  straight  forward,  however,  to  compare  two  sparse 
spectra  because  1)  the  number  of  relaxations  can  be  different  and  2)  small  misalignments 
of  two  relaxation  frequencies  often  occur.  For  these  two  reasons,  the  standard  Euclidean 
distance  is  not  an  ideal  distance  measure  between  two  DSRFs.  For  example,  in  Fig.  10(a) 
the  estimate  and  the  correct  solution  are  very  similar  but  with  a  large  Euclidean  distance 
of  0.71  units  of  c  because  of  one  small  misalignment  in  the  (  axis.  This  issue  is  further 
complicated  when  the  number  of  relaxations  in  two  spectra  are  not  the  same.  If  two  DSRFs 
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S  and  S  have  the  same  model  order,  the  Euclidean  distance  between  the  two  is 

—2 

(19) 

The  Earth  Mover’s  Distance  (EMD)  is  a  measure  that  takes  into  account  misalignments 
and  mismatching  model  orders  [35,43].  The  EMD  measures  how  much  “work”  it  takes  to 
morph  one  spectrum  into  another.  The  metaphor  is  that  one  spectrum  represents  piles  of 
earth  with  volume  c*,  located  at  the  associated  Ck-  The  other  spectrum  represents  holes  in 
the  ground  with  capacity  located  at  Ck-  The  EMD  is  proportional  to  the  least  amount 
of  work  needed  to  move  as  much  earth  as  possible  into  the  holes.  Work  is  defined  to  be  the 
amount  of  earth  moved  times  the  distance  traveled.  Intuitively,  the  distance  would  simply 
be  the  difference  between  Ck  and  Ck  in  log-space.  This  is  true  when  Ck  and  Ck  have  the  same 
sign.  To  accommodate  opposite  signs,  the  distance  function  dij  between  two  relaxations 
(Cj,  Cj)  and  (Cj,  dj )  is  defined  to  be 


K 


^  '  (Cfc  Ck) 


k=l 


dij 


I  log  0  —  log  0 1  ,  CiCj  >  0 

1  +  I  logCi  -  log  01  ,  CiCj  <  0 


which  penalizes  relaxations  with  different  signs  by  adding  a  constant  to  the  distance.  This 
constant  can  be  replaced  by  other  quantities,  such  as  exp(— a|log^  —  log  C?  I ) ,  where  a 
determines  the  decaying  rate.  Spectra  are  made  nonnegative  and  normalized  prior  to  the 
EMD  computation.  The  EMD  is  measured  in  decades  because  it  is  examined  in  log-C  space. 

The  EMD  compares  two  spectra  as  a  whole,  so  the  effect  of  very  small  amplitude  re¬ 
laxations  is  tiny  in  the  EMD,  and  neglecting  these  small  components  amounts  to  assuming 
they  are  near  the  noise  level  of  the  measured  frequency  response.  See  Appendix  C  for  details 
about  the  EMD.  For  simplicity,  it  is  sufficient  to  say  that  EMD  quantifies  the  similarity 
between  two  DSRFs. 

Some  examples  of  the  EMD  are  shown  in  Fig.  10.  It  is  observed  that  the  EMD  reflects 
the  similarity  between  two  DSRFs.  On  the  other  hand,  the  £ 2  Euclidean  distance  returns  a 
similar  number  regardless  of  how  similar  two  DSRFs  are  to  each  other. 
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Figure  10:  Examples  of  the  £2  Euclidean  distance  and  the  EMD  between  two  DSRFs, 
labeled  by  the  diamonds  and  the  squares,  (a)  The  two  spectra  are  very  similar  but  the  £2 
distance  is  large  while  the  EMD  is  small,  (b)  The  spectra  are  less  similar  compared  to  (a), 
which  is  reflected  in  the  higher  EMD,  but  the  £2  distance  remains  the  same,  (c)  Spectra 
with  different  model  orders,  (d)  Two  very  different  spectra. 


2.2.2  Synthetic  Six-relaxation  DSRF 

The  proposed  sparsity-promoting  estimation  method  is  tested  on  a  six-relaxation  DSRF 
synthesized  at  70 dB  SNR  with  additive  white  Gaussian  noise  (Fig.  11).  This  is  a  case  that 
cannot  be  handled  by  traditional  nonlinear  parameter  optimization  [30] ,  or  the  nonnegative 
linear  method  [33].  The  actual  DSRF  and  the  frequency  response,  along  with  their  fits,  are 
shown  in  Fig.  11.  All  six  relaxation  frequencies  are  recovered  by  solving  (18)  with  p  =  0.5. 
The  estimation  is  nearly  perfect,  because  the  estimated  model  parameters  are  real,  and  the 
deviation  from  the  true  answer  is  small;  the  EMD  is  0.09  decade,  for  p  =  0.5. 

When  this  spectrum  is  estimated  with  p  =  1  using  ll_ls  many  extra  relaxations  are 
introduced  by  the  fitting  process;  for  p  =  1,  the  EMD  is  0.10  decades.  Real  targets  are  not 
likely  to  have  a  spectrum  with  so  many  small  relaxations  around  a  strong  relaxation.  In 
fact,  Baum  argues  that  physical  relaxation  frequencies  are  discrete  [12],  However,  the  several 
small  relaxations  when  using  p  =  1  give  a  denser  distribution  of  relaxations,  resembling  a 
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continuous  spectrum.  In  this  sense,  p  <  1  gives  a  sparser  solution  that  more  accurately 
resembles  a  physical  spectrum  even  though  this  may  not  always  be  reflected  in  the  EMD 
measure.  The  fits  of  the  frequency  response  is  also  displayed  in  Fig.  11(a),  and  it  is  seen 
that  both  p  =  0.5  and  p  =  1  produce  a  good  fit  to  the  frequency  response. 
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Figure  11:  Estimation  of  a  synthetic  six-relaxation  DSRF  with  p  =  0.5  and  p  =  1.  (a)  The 
synthesized  frequency  response  and  two  model  fits,  (b)  The  actual  DSRF  and  its  estimates. 

2.2.3  Signal-to-Noise  Ratio 

To  see  how  the  proposed  method  performs  in  noise,  a  Monte  Carlo  simulation  of  the  goodness 
of  estimation  versus  SNR  is  performed.  For  each  SNR,  100  samples  are  generated,  and 
each  sample  has  a  DSRF  with  four  relaxations.  The  relaxation  frequencies  are  uniformly 
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distributed  at  random  over  [Cmin,  Cm  ax] ,  and  the  spectral  amplitudes  are  generated  similarly 
over  [—1,  1].  For  each  DSRF  estimate,  the  goodness  of  estimation  is  measured  by  the  EMD 
between  the  estimate  and  truth. 

The  simulation  result,  in  Fig.  12,  shows  the  robustness  of  the  estimation  method  at 
different  SNRs.  The  EMD  between  the  estimate  and  the  truth  increases  as  the  SNR  de¬ 
creases.  This  suggests  that  the  proposed  method  is  usable  in  a  range  of  SNR  where  the 
EMD  is  below  some  threshold.  Visually,  spectra  with  an  EMD  below  0.1  are  considered 
similar,  those  with  an  EMD  above  0.2  exhibit  noticeable  differences,  but  may  still  resemble 
each  other.  In  our  laboratory  measurements,  a  typical  SNR  for  loop  targets  is  70  dB. 


Figure  12:  Monte  Carlo  simulation  of  goodness  of  estimation  (EMD)  vs.  SNR  performed 
using  a  four-relaxation  DSRF.  The  error  bars  indicate  the  range  of  EMD  between  the  10th 
and  90th  percentiles.  Sample  size  is  100  per  SNR. 

2. 3  Laboratory  Data 

To  verify  the  functionality  of  the  algorithm,  the  proposed  method  is  applied  to  the  lab¬ 
oratory  measurement  of  a  target  where  the  theoretical  DSRF  is  known  [44],  The  target 
consists  of  three  mutually-orthogonal  copper  loops.  The  loop  diameters  and  thicknesses  are 
3/20,  4/30,  and  5/36,  respectively  in  crn/AWG1.  The  orientation  and  position  relative  to 
the  EMI  sensor  are  chosen  to  demonstrate  the  positive  and  negative  spectral  amplitudes. 
The  measured  frequency  response  of  this  configuration  is  shown  in  Fig.  13(a),  and  the  SNR 
is  about  38 dB.  The  estimated  DSRF  is  displayed  in  Fig.  13(b).  The  estimate  and  theory 
agree  well,  and  the  EMD  between  the  theoretical  and  estimated  DSRF  is  0.10  decades. 

1  American  Wire  Gauge 
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Figure  13:  (a)  Measured  frequency  response  of  three  mutually  orthogonal  loops,  its  fit, 
and  the  theoretical  response,  (b)  Theoretical  and  estimated  DSRFs  of  the  response  in  (a). 

Next,  changes  in  the  DSRF  as  the  target  moves  relative  to  the  EMI  sensor  are  studied. 
The  same  target  configured  at  a  fixed  orientation  is  displaced  at  different  positions  along 
a  horizontal  axis,  denoted  x.  The  vertical  distance  between  the  target  and  sensor  is  6  cm. 
The  EMI  sensor  is  located  at  x  =  0.  Samples  of  the  measured  target  responses  are  shown 
in  Fig.  14(a);  their  corresponding  spectra  are  in  Fig.  14(b).  Theoretical  results  are  also 
shown.  Overall,  the  theory  and  measurement  agree.  The  disagreement  at  x  =  —0.5  may  be 
because  of  approximations  in  the  model  and/or  inaccuracies  in  the  positions  measured  in 
the  experiment. 

As  expected  from  the  theory,  while  the  frequency  response  changes  dramatically  as  the 
target  moves  along  the  x  axis,  the  corresponding  change  in  the  spectral  domain  only  occurs 
in  the  spectral  amplitudes.  The  three  dominant  relaxation  frequencies  remain  unchanged. 
The  proposed  method  successfully  estimates  the  spectra  that  agree  with  this  phenomenon. 
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Figure  14:  The  plots  share  the  same  annotation  as  Fig.  13.  (a)  Frequency  responses 
of  the  three  mutually  orthogonal  copper  loops  at  different  locations,  (b)  Theoretical  and 
estimated  DSRFs  of  the  corresponding  responses  in  (a).  The  SNR  is  measured  in  dB,  x 
positions  in  cm,  and  EMD  in  decades. 


All  three  relaxation  frequencies  are  consistently  estimated.  The  extra  relaxations  all  have 
small  amplitudes  that  can  be  safely  ignored.  This  invariant  property  of  the  relaxation 
frequencies  makes  the  DSRF  valuable  especially  for  target  discrimination. 
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2.^  Field  Data 


The  proposed  method  is  applied  to  a  data  set,  acquired  during  a  field  test,  which  contains 
62  types  of  targets,  including  26  types  of  landmines  as  well  as  various  types  of  metallic  and 
nonmetallic  clutter.  The  testing  field  is  divided  into  11  lanes,  each  lane  containing  20  grid 
cells,  so  there  are  220  grid  cells  total.  About  145  measurements  are  collected  per  grid. 

The  primary  objective  of  estimating  the  DSRF  from  field  measurements  is  to  verify  the 
functionality  of  the  proposed  method  for  realistic  field  data,  and  to  study  the  consistency 
of  DSRFs  from  targets  of  the  same  type.  Because  targets  of  the  same  type  should  have  the 
same  physical  composition,  their  DSRFs  are  expected  to  be  the  same.  It  is  demonstrated  in 
the  following  section  that  targets  from  the  same  type  do  have  very  similar  estimated  DSRFs. 
Targets  of  different  types  are  also  observed  to  have  very  different  DSRFs,  in  general. 

The  field  measurements  are  collected  by  an  EMI  measurement  cart,  shown  in  Fig.  1 
in  Chapter  1.  The  EMI  acquisition  hardware  used  has  the  same  specification  as  described 
in  Section  2.2,  and  details  are  found  in  [8].  As  the  cart  is  pushed  over  a  target,  the  fre¬ 
quency  response  of  the  target  is  collected  sequentially,  as  illustrated  in  Fig.  2.  A  segment 
of  collected  responses  (magnitude  only)  at  21  frequencies  is  shown  in  Fig.  15.  Although 
multiple  measurements  are  available  per  target,  in  this  chapter  only  one  measurement  is 
considered  per  target.  The  one  measurement  selected  has  the  largest  response  magnitude 
in  the  grid  cell.  In  Fig.  15,  this  corresponds  to  the  peak  of  the  strongest  lobe  in  a  grid  cell, 
which  usually  occurs  when  the  EMI  sensor  is  directly  above  or  nearly  above  the  target.  The 
multiple  lobes  seen  in  Fig.  15  are  a  result  of  filtering  the  measurements  [8].  For  field  mea¬ 
surements,  the  SNR  is  estimated  by  the  signal-to-background  ratio,  where  the  background 
is  considered  to  be  primarily  due  to  the  soil.  The  field  measurements  used  in  this  work  have 
a  background  signal  power  of  about  — 130  dB,  as  shown  in  grid  70  in  Fig.  15. 

In  field  measurements,  the  soil  response  is  also  measured  in  addition  to  the  target 
response  and  noise.  In  Section  2.4.2,  the  frequency  response  of  soil  and  its  effect  on  the 
DSRF  are  examined.  An  “augmented  dictionary”  that  includes  the  DSRF  model  and  a  soil 
model  is  proposed  to  simultaneously  model  the  target  response  and  the  soil  response.  From 
preliminary  results,  it  is  observed  that  including  the  soil  model  in  the  DSRF  estimation 
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process  produces  estimates  very  similar  to  those  obtained  when  ignoring  the  soil  response. 
This  observation  suggests  that  it  is  safe  to  ignore  the  soil  response  when  estimating  the 
DSRF  of  targets  from  held  measurements.  However,  more  simulations  and  tests  should  be 
done  to  study  the  effect  of  soil  on  target  responses.  The  results  presented  in  the  following 
section  are  obtained  ignoring  the  effect  of  soil. 


Figure  15:  A  short  segment  of  the  held  measurements  acquired  by  the  EMI  cart.  Shown 
are  the  magnitude  of  the  frequency  responses;  each  curve  represents  one  frequency.  The 
letters  denote  the  type  of  target  buried  in  the  grid  cell. 

2.4.1  Estimating  the  DSRF  of  Various  Targets 

Several  types  of  landmines  are  examined  to  study  the  performance  of  the  proposed  method 
when  applied  to  held  measurements  and  to  investigate  the  consistency  of  DSRF  from  land¬ 
mines  of  the  same  type.  Each  type  of  landmine  has  multiple  instances  buried  at  different 
depths  in  different  grid  cells.  The  buried  depth  affects  the  strength  of  the  measured  re¬ 
sponse.  One  sample  of  the  measurements  is  used  per  instance,  and  the  DSRF  of  each 
sample  is  estimated  and  then  plotted  together  with  others  of  the  same  type.  The  spectral 
amplitudes  are  represented  by  the  color  intensity.  The  DC  magnetization  of  targets  is  also 
examined  at  the  end  of  this  section. 

Eight  Type-A  mines,  a  low-metal  content,  nonmagnetic,  and  moderate  EMI  response 
antipersonnel  mine,  are  examined.  The  SNR  ranges  from  about  45  to  60  dB.  The  fre¬ 
quency  response  of  the  eight  instances  are  shown  in  Fig.  16(a)  and  the  estimated  DSRFs  in 
Fig.  16(b).  All  eight  Type-A  mines  exhibit  consistency  in  the  relaxation  frequencies  and  the 
spectral  amplitudes.  The  DSRFs  are  very  similar,  and  the  average  EMD  between  pairs  of 
mines  is  0.052  decades.  The  estimated  DSRFs  also  generate  good  curve  fits  in  the  frequency 
domain,  shown  in  Fig.  16(a). 
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Eight  Type-B  mines,  a  medium-metal  content,  magnetic,  strong  EMI  response  antiper¬ 
sonnel  mine,  are  examined.  The  SNR  ranges  from  about  55  to  70  dB.  The  frequency  re¬ 
sponses  are  shown  in  Fig.  16(c)  and  the  estimated  DSRFs  in  Fig.  16(d).  The  samples  exhibit 
consistency  in  the  relaxation  frequencies  and  the  spectral  amplitudes.  Samples  #1  and  #3 
differ  from  other  samples  slightly  near  3.7  decades  where  two  relaxations  are  combined  into 
one,  but  the  difference  is  small.  The  average  EMD  between  pairs  of  mines  is  0.1  decades. 
The  estimated  DSRFs  generate  good  curve  fits  in  the  frequency  domain. 

To  demonstrate  the  advantage  of  the  sparsity-promoting  DSRF  estimation  method  over 
the  NNLSQ,  four  samples  of  the  Type-K  mine,  a  medium-metal  content,  magnetic,  moderate 
EMI  response  antipersonnel  mine,  are  fitted  using  both  methods.  The  SNR  ranges  from 
about  52  to  60  dB.  The  frequency  responses  and  the  estimated  DSRFs  are  shown  in  Fig.  17. 
From  the  frequency  response  plots,  it  is  observed  that  sample  #4  has  a  relaxation  with 
a  negative  Ck  because  the  real  part  of  the  response  is  non-monotonic.  This  property  was 
discovered  during  an  early  phase  of  this  research  [34]  and  is  used  as  the  basis  for  using 
the  NNLSQ  because  most  filtered  field  target  responses  exhibit  a  monotonic  real  part  [33]. 
However,  sample  #4  does  not  follow  this  trend  and  has  a  negative  component.  The  NNLSQ 
method  returns  a  DSRF  estimate  for  #4  that  is  quite  different  from  other  samples,  and  the 
fit  in  the  frequency  response  is  also  not  very  good.  The  low  frequency  part,  where  the  real 
part  is  non-monotonic,  is  not  fitted  by  NNLSQ.  On  the  other  hand,  the  sparsity-promoting 
method  property  recovers  the  negative  relaxations  in  $4  and  the  estimated  DSRF  agrees 
with  other  samples. 

More  examples  of  estimated  DSRFs  from  field  targets  are  shown  in  Fig.  18  and  Fig.  19. 
Here  is  a  list  of  description  for  the  example  targets: 

•  Type-H  and  -I:  low-metal  content,  magnetic,  moderate  response  antitank  mines. 

•  Type-V  and  -W:  low-metal  content,  magnetic,  moderate  response  antipersonnel  mines. 

•  Type-C:  a  low-metal  content,  magnetic,  weak  response  antipersonnel  mine. 

•  Type-D:  a  low-metal  content,  slightly  magnetic,  strong  response  antipersonnel  mine. 

•  Type-E:  a  low-metal  content,  magnetic,  weak  response  antipersonnel  mine. 

•  Type-L:  a  medium-metal  content,  magnetic,  strong  response  antipersonnel  mine. 
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Consistent  DSRFs  within  a  target  type  are  generally  observed.  Deviations  do  occur,  but 
mostly  in  relaxations  with  a  small  spectral  amplitude,  e.g.,  Type-E  mines  in  Fig.  19(c).  A 
small  spectral  amplitude  contributes  only  a  small  amount  to  the  EMD.  Another  observation 
is  the  combining  or  splitting  of  relaxations,  such  as  the  Type-W  mines  in  Fig.  18(d).  Often 
the  combining  of  relaxations  happens  at  lower  SNRs,  where  there  is  not  enough  resolution 
to  separate  the  two  relaxations.  This  is  reflected  in  the  A  selection,  where  a  lower  SNR 
gives  a  higher  A,  encouraging  a  sparser  solution,  as  discussed  in  Section  2.5.  Nevertheless, 
the  impact  of  combining  and  splitting  is  inherently  reduced  in  the  EMD.  Targets  with  a 
weak  EMI  response,  such  as  Type-C  mines  in  Fig.  19(a),  are  also  observed  to  have  more 
variations  in  the  estimated  DSRF  and  higher  EMDs. 

On  the  other  hand,  targets  of  different  types  tend  to  have  very  distinct  DSRFs.  For  ex¬ 
ample,  when  comparing  Type-A  and  Type-B  mines,  Type-A  mines  only  have  two  relaxations 
and  Type-B  mines  have  six  relaxations.  Even  if  the  number  of  relaxations  are  the  same,  the 
spectral  amplitudes  can  be  different.  For  example,  Type-H  and  Type-V  mines  both  have 
four  relaxations  but  Type-H  has  its  stronger  relaxations  around  6  decades,  whereas  Type-V 
has  its  strongest  relaxation  near  5  decades.  There  are  mines  from  the  same  family  that  do 
have  similar  DSRFs,  such  as  Type-V  and  Type-W  mines.  The  characteristics  of  each  target 
type  observed  in  the  DSRF  suggest  a  DSRF-based  target  discrimination,  which  is  discussed 
in  Chapter  4. 

Estimates  of  H{ 0)  =  Ylk= o  f°r  several  mine  types  are  shown  normalized  in  Fig.  20. 
This  constant  is  due  to  the  DC  magnetization  of  the  target.  For  example,  Type-A  mines  are 
nonmagnetic  and  the  estimated  H( 0)  is  close  to  zero,  as  shown  in  Fig.  20(a).  In  contrast, 
Type-V  and  -W  mines  are  strongly  magnetic  and  have  a  relatively  large  H( 0).  Type-I 
and  -B  are  magnetic  but  have  a  weaker  magnetic  response  than  Type-V  and  -W.  Shown 
in  Fig.  20(b)  are  some  cases  with  poor  H( 0)  estimates.  For  example,  sample  #3  of  the 
Type-H  mine  has  a  weak  target  response  and  an  out-of-band  relaxation  that  is  not  fitted 
well,  making  cfc  a  poor  estimate  for  0).  Sample  #5  of  the  Type-C  mine  also  has  a 

very  weak  target  response  at  18  dB  SNR,  and  H( 0)  is  likely  corrupted  by  the  soil. 


31 


Figure  16:  Frequency  response  and  estimated  DSRF  of  Type-A  (a,b)  and  Type-B  (c,d)  mines  measured  in  the  field,  (a)  The  frequency 
response  of  eight  Type-A  mines  (solid  lines)  and  their  fits  (square  markers),  (b)  Estimated  DSRFs  from  the  eight  frequency  responses, 
(c)  The  frequency  response  of  eight  Type-B  mines  and  their  fits,  (d)  Estimated  DSRFs. 


imag  x  imag 


Figure  17:  Comparison  of  the  proposed  sparsity  method  (a,b)  and  the  NNLSQ  (c,d).  The  measured  frequency  responses  (solid  lines) 
and  their  fits  (square  markers)  are  shown  in  (a)  and  (c),  and  the  estimated  DSRF  in  (b)  and  (d).  Colors  red  and  blue  are  used  to 
represent  positive  and  negative  spectral  amplitudes,  respectively. 
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(c) 


(d) 


Figure  19:  Examples  of  estimated  DSRFs  from  field  targets.  The  average  EMD  listed  below  are  in  decades  and  SNR  in  dB.  (a)  Eight 
Type-C  mines;  EMD=0.2,  SNR=26.  (b)  Eight  Type-D  mines;  EMD=0.085,  SNR=82.  (c)  Eight  Type-E  mines;  EMD=0.093,  SNR=26. 
(d)  Four  Type-L  mines;  EMD=0.043,  SNR=85. 


Figure  20:  Normalized  estimated  H( 0)  for  several  mine  types.  The  normalization  factor 
is  ^2 Cfc.  (a)  Good  and  consistent  H( 0)  estimates,  (b)  Mine  types  with  instances  of  poor 
H( 0)  estimates  due  to  low  SNR  and  out-of-band  relaxations. 


2.4.2  Soil  Response  and  Soil  Model 

To  study  the  frequency  response  of  soil  and  its  effect  on  the  DSRF,  soil  responses  from  the 
field  are  collected  and  examined.  These  measured  frequency  responses  are  collected  at  grid 
cells  reported  to  have  no  objects  buried,  such  as  grid  70  in  Fig.  15.  The  frequency  response 
of  a  subset  of  these  soil  responses  are  plotted  in  Fig.  21.  It  is  observed  that  the  frequency 
dependence  of  the  soil  responses  share  a  similar  trend.  The  real  part  has  a  linear  trend  with 
respect  to  the  log- frequency,  and  the  imaginary  part  tends  to  be  a  constant  [8,45]. 

From  these  observations,  the  following  soil  model  is  suggested: 

HG(u)  =  'y(]n-+jf),  (20) 

V  w0  2  J 

where  7  and  ujq  are  model  parameters.  More  details  about  this  model  can  be  found  in  [46]. 
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Figure  21:  Instances  of  measured  EMI  response  of  soil,  (a)  The  real  part  has  a  linear 
trend  with  respect  to  the  log-frequency,  and  (b)  the  imaginary  part  tends  to  be  a  constant. 


To  fit  a  measured  response  to  the  soil  model,  (20)  is  rewritten  as 

Hg(u)  =Pi  +P2  (hiw  +  ,  (21) 

where  pi  =  —7  lncco  and  p2  =  7  are  model  parameters.  For  a  soil  response  measured  at 
N  frequencies, 

&G  =  Gp  +  noise  ,  (22) 

where 

1  Inter  +  j  7t/2 

1  lncj2  +  J7r/2  Pi 

G  =  and  p  = 

:  :  P2 

1  lno;jv+J7r/2 

To  consider  the  target  response  in  the  presence  of  soil,  as  described  in  Section  1.1,  the 
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total  response  &Totai  can  be  modeled  as 


&Totai  =  b  +  bG  +  noise,  (23) 

where  b  is  the  target  response.  Combining  the  DSRF  model  (16)  and  the  soil  model  (22), 
the  model  (23)  can  be  written  as  a  set  of  linear  equations: 

&Totai  =  Ax  +  Gp  +  e,  (24) 

where  e  accounts  for  the  modeling  error  as  well  as  the  measurement  noise.  The  two  matrices 
A  and  G  can  then  be  consolidated  into  an  augmented  coefficient  matrix  <1>.  That  is, 

(25) 


(26) 


in  the 

augmented  form,  represented  by  the  first  entry  of  0. 

Using  the  augmented  model  (25),  the  target  response  and  the  soil  response  can  be 
modeled  simultaneously.  The  approach  to  solving  (25)  is  identical  to  that  of  solving  the 
DSRF  modeling  problem  (18),  that  is,  minimize  the  squared  error  ||6xotal  —  ^0111  subject 
to  a  sparsity  regularization  ||0||p.  The  sparsity  regularization  works  for  the  augmented 
problem  because  the  sparsity  of  the  solution  increases  only  by  one,  to  account  for  the  second 
column  of  <F,  which  corresponds  to  the  soil  model.  Therefore,  the  solution  0  should  still  be 
sparse.  In  short,  the  target  response  and  the  soil  response  can  be  modeled  simultaneously 
by  solving  the  sparsity-regularized  least-squares  (18)  with  the  dictionary  A  replaced  by  the 
augmented  matrix  $. 


-'Total 


=  $0  +  e, 


where 


$  = 


1  lnwi  +  j7r/2 
1  lnw2+Jvr/2 

1  lnLUN+jn/2 


1-JW1/C2 

1— JWi/Cm 

i+i^i/6 

1+JW1/C2 

i+jWCm 

1 -1^2 /Cl 

1— JW2/C2 

1—  3^2 /Cm 

l+jW2/Cl 

1+/W2/C2 

l+jlV2/CM 

1— JWjv/ Cl 

1—  j^>N  K2 

1— /^iv/Cm 

i+hWCi 

1+JWjv  /  C2 

I+I^iv/Cm  _ 

and  0  = 


do  +  Pi 
P2 
ci 
C2 

cm 


The  shift  coefficient  in  the  DSRF  model  cq  and  the  soil  model  p\  are  combined 
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The  effect  of  the  soil  on  the  field  target  response  is  considered  here.  This  is  done  by 
comparing  DSRF  estimates  obtained  with  the  soil  model  via  the  augmented  form  (25)  and 
without  the  soil  model  (16).  Several  landmines  are  examined,  and  it  is  found  that  the  two 
approaches  give  very  similar  and  sometimes  identical  estimates.  Four  examples  are  given 
in  Fig.  22  and  Fig.  23. 

In  Fig.  22,  eight  Type-A  and  eight  Type-B  mines  are  examined.  For  Type-A  mines, 
including  the  soil  model  does  not  produce  different  DSRF  estimates.  For  Type-B  mines, 
including  the  soil  model  produces  DSRF  estimates  that  are  very  similar  to  those  obtained 
when  not  using  the  soil  model.  Only  slight  shifts  in  the  relaxation  frequencies  are  observed. 
The  averaged  EMD  between  the  two  sets  is  0.023  decades,  which  is  computed  by  averaging 
over  the  instances  the  EMD  between  the  DSRFs  obtained  with  and  without  the  soil  model. 

In  Fig.  23,  eight  Type-C  and  eight  Type-D  mines  are  examined.  For  Type-C  mines, 
including  the  soil  model  produces  almost  identical  DSRF  estimates.  The  average  EMD 
between  the  two  sets  is  0.0082  decades.  For  Type-D  mines,  including  the  soil  model  produces 
DSRF  estimates  that  are  very  similar  to  those  obtained  not  using  the  soil  model.  A  slight 
difference  is  observed  in  relaxation  frequencies  near  3.5  and  6  decades.  The  averaged  EMD 
between  the  two  sets  is  0.027  decades. 

The  four  types  of  mines  shown  in  Fig.  22  and  23  cover  a  wide  range  of  SNR,  from 
18  dB  to  71  dB,  and  the  consistency  between  the  DSRFs  estimated  with  or  without  the 
soil  model  is  observed  in  all  cases.  The  results  show  that  estimating  the  DSRF  without 
the  soil  model  (16)  produces  a  DSRF  that  is  very  similar  to  that  obtained  with  the  soil 
model  (25).  An  explanation  for  this  observation  is  that  the  soil  model  does  not  fit  well  to 
a  sparse  DSRF  model,  so  it  is,  in  a  sense,  orthogonal  to  the  DSRF  terms.  Therefore,  not 
including  the  soil  term  in  the  modeling  does  not  result  in  the  soil  response  “leaking”  into 
the  DSRF  terms.  Another  reason  is  simply  that  the  soil  response  is  negligible  when  the 
target  response  is  much  stronger  than  the  soil  response.  In  either  case,  the  results  suggest 
that  it  is  safe  to  ignore  the  soil  in  estimating  the  DSRF  from  field  measurements.  However, 
further  simulations  should  be  performed  to  study  and  quantify  the  range  of  SNR  over  which 
the  soil  can  be  ignored. 
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Sample  Sample 


Figure  23:  Comparisons  of  the  estimated  DSRF  when  including  the  soil  model  (a,b)  and  not  including  the  soil  model  (c,d).  Shown  in 
(a)  and  (c)  are  Type-C  mines,  and  the  average  EMD  is  0.082  decades  between  the  two  sets.  Shown  in  (b)  and  (d)  are  Type-D  mines, 
and  the  average  EMD  is  0.027  decades  between  the  two  sets. 


2.5  Choosing  the  Regularization  Parameter  A 


In  this  section,  it  is  first  examined  the  behavior  of  the  proposed  method  in  relation  to  the 


regularization  parameter  A,  and  then  a  simple  A  selection  rule  is  proposed  exploiting  the 
observed  properties  of  A.  All  discussions  and  figures  presented  here  assume  p  =  0.5  unless 
otherwise  specified. 

To  understand  how  the  goodness  of  fit  changes  with  A  and  SNR,  a  cross-validation¬ 
like  simulation  is  conducted.  First,  a  collection  of  synthetic  spectra  are  built  with  different 
model  orders  and  a  variety  of  distributions  of  relaxations.  For  each  spectrum  at  a  fixed  SNR, 
the  spectrum  is  estimated  100  times  for  each  A  within  a  range,  and  the  average  goodness 
of  fit,  measured  by  the  EMD  between  the  available  truth  and  the  estimate,  is  recorded. 
This  is  done  for  a  range  of  SNRs.  The  simulation  result  for  a  four-relaxation  spectrum,  as 
an  example,  is  shown  in  Fig.  24.  Not  only  is  the  EMD  surface  well-behaved  (i.e. ,  smooth) 
with  respect  to  the  SNR  and  A,  but  more  importantly  the  surface  itself  is  convex-shaped. 
Thus,  at  each  SNR,  the  minimum  EMD  is  achievable  by  a  unique  A.  The  wide  valley  of 
the  surface  also  shows  that  the  goodness  of  fit  is  not  very  sensitive  near  the  optimal  A  that 
gives  the  minimum  EMD  per  SNR. 


x 


SNR(dB) 


Figure  24:  Monte  Carlo  simulation  of  the  goodness  of  estimation  (EMD)  of  a  four- 
relaxation  spectrum  at  different  SNR’s  and  A’s. 
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Simulations  of  spectra  for  other  model  orders  and  distributions  also  exhibit  the  same 
property  (Fig.  25).  Moreover,  the  valleys  of  the  EMD  surface  all  occur  in  nearly  the  same 
SNR-A  region.  In  other  words,  the  A  that  produces  the  minimum  EMD  at  a  given  SNR 
is  quasi-independent  of  the  model  order.  Figure  26  shows  the  averaged  EMD  of  different 
model  orders  in  Fig.  25.  The  resulting  surface  still  exhibits  the  properties  described  above. 
This  allows  us  to  pick  a  near-optimal  A  based  solely  on  the  SNR. 


Figure  25:  Same  simulation  as  in  Fig.  24,  but  for  spectra  of  different  model  orders  (1  to 
6).  Each  spectrum  constitutes  one  surface  in  the  figure.  All  surfaces  have  their  minimum 
in  the  same  SNR-A  region. 

In  Fig.  26,  the  optimal  A  at  each  SNR  which  is  also  plotted.  Using  the  wide-valley 
property,  the  near-minimum  EMD  can  be  achieved  by  choosing  A’s  that  are  near  the  optimal 
A  in  the  valley.  Here  the  optimal  A  is  approximated  with  a  semilog  function  of  SNR.  This 
is  done  by  fitting  the  optimal  log-A  curve  with  a  linear  function.  Weights  may  be  added  to 
promote  certain  SNR’s  that  are  more  important.  For  our  problem  setup,  the  A  is  chosen  by 
(also  shown  in  Fig.  26) 

log  A  =  -0.05 -SNR -2.2  (27) 

In  practice,  this  log-A  selection  rule  that  is  linear  in  SNR  allows  the  regularization  parameter 
to  be  determined  with  negligible  computation  time.  When  processing  the  laboratory  data, 
(27)  is  used  along  with  an  estimate  of  the  SNR  to  determine  A  for  use  in  Algorithm  1. 
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The  same  empirical  method  can  be  repeated  for  other  p’s,  and  the  result  is  also  a  linear 
relationship  between  log  A  and  SNR. 

Projection  of  the  optimal  A.  path 


Figure  26:  Average  of  EMD  surfaces  in  Fig.  25.  The  curve  with  asterisk  markers  traces 
out  the  optimal  A’s.  The  line  with  square  markers  is  the  approximated  optimal  A  curve 
used  to  select  A  in  practical  estimation. 

Figure  27  compares  the  goodness  of  estimation  of  a  four-relaxation  target  using  the  linear 
log-A  selection  rule  and  the  optimal  A  which  uses  the  true  spectrum  that  is  not  available  in 
practice.  There  is  a  slight  increase  in  the  EMD  when  the  linear  log-A  selection  rule  is  used, 
which  is  reasonable  and  expected.  The  increase  is  acceptable  and  hence  the  linear  log-A 
selection  rule  is  an  appropriate  A  selector. 

Also  shown  in  Fig.  27  are  the  performances  of  other  p  values.  It  is  seen  that  p  <  1 
gives  more  accurate  results  than  p  =  1  when  the  optimal  A  is  used,  but  this  advantage  is 
significantly  diminished  when  the  linear  log-A  selection  rule  is  used.  While  this  lessens  the 
advantages  of  using  p  <  1,  it  is  emphasized  that  p  =  1  tends  to  give  estimates  with  many 
relaxation  frequencies  while  p  <  1  gives  sparser  spectra  which  are  more  physically  accurate 
(see  Section  2.2.2).  It  is  possible  that  both  the  accuracy  and  the  sparsity  advantages  for 
p  <  1  could  be  obtained  with  a  better  A  selection  rule.  Lastly,  since  the  performance  of  a 
certain  p  value  is  dependent  on  the  A  selection  rule  used,  different  optimum  p  values  would 
be  determined  if  the  A  selection  rule  is  changed. 
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Figure  27:  Goodness  of  fit  using  the  linear  log-A  selection  rule  for  several  p' s.  The  true 
spectrum  is  the  same  as  in  Fig.  24.  Dash-dot  curves  denote  the  optimal  A,  solid  lines  denote 
the  linear  log-A  selection  rule. 
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CHAPTER  III 


SPECTRUM  ESTIMATION  FROM  MULTIPLE  MEASUREMENTS 


In  the  process  of  acquiring  the  EMI  response  of  a  target,  often  the  target  response  is 
measured  at  different  target-to-sensor  orientations  and  locations.  For  example,  this  happens 
when  an  EMI  sensor  collects  data  while  moving  over  a  target,  resulting  measurement  views 
of  the  target  from  different  angles,  as  illustrated  in  Fig.  2  in  Chapter  1.  It  seems  reasonable 
that  more  accurate  estimates  of  the  DSRF  can  be  obtained  with  more  views  or  snapshots 
of  the  target,  for  two  reasons:  1)  the  relaxation  frequencies  are  constant  in  different  views 
and  2)  a  relaxation  that  vanishes  (c*,  =  0)  in  one  view  can  appear  in  a  different  view,  as 
illustrated  in  Fig.  13  in  Chapter  2. 

To  utilize  the  multiple  measurements  often  available  for  a  target,  a  multiple- measurement 
vector  (MMV)  DSRF  estimation  method  is  developed  in  this  chapter.  The  MMV  method 
exploits  the  property  of  orientation  and  position  invariance  of  the  relaxation  frequencies, 
and  obtains  more  accurate  estimates  by  encouraging  different  measurements  (views  of  a 
target)  that  share  a  common  set  of  relaxation  frequencies.  The  proposed  MMV  method  is 
applied  to  synthetic,  laboratory,  and  field  data,  and  the  performance  is  demonstrated  to  be 
robust  and  better  than  the  SMV  approach. 

The  MMV  DSRF  estimation  is  a  generalization  of  the  single-measurement  vector  (SMV) 
approach  (Chapter  2),  but  the  MMV  generalization  does  more  because  the  SMV  does  not 
exploit  the  invariance  of  the  relaxation  frequencies.  The  MMV  has  an  advantage  over  the 
SMV  by  taking  into  account  this  physical  invariance  property. 

A  “row-sparsity  measure”  for  matrices  is  introduced  to  encourage  DSRF  estimates  for 
the  same  target  to  share  the  same  relaxation  frequencies.  This  measure  formulates  the 
MMV  problem  in  a  form  similar  to  the  SMV,  allowing  SMV  techniques  to  be  extended  to 
the  MMV  case  in  a  straightforward  manner. 
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The  MMV  method  follows  the  same  framework  as  the  SMV:  1)  sample  the  relaxation  fre¬ 
quency  domain,  2)  construct  an  overcomplete  dictionary,  3)  perform  a  sparsity-promoting 
optimization,  and  4)  interpolate.  Because  steps  1)  and  2)  are  identical  to  the  SMV  ap¬ 
proach,  the  discussion  here  is  focused  on  the  MMV  sparsity-promoting  optimization  and 
interpolation. 

3.1  Method  Formulation 


where  B  G  CNxL  has  columns  of  b±, . . . ,  b i,  and  similarly  for  X  G  ]g>(A/+i)xL  and  E  G 
CNxL.  If  the  L  measurements  come  from  the  same  target,  then  x\, . . .  ,xl  share  the  same 
location  of  nonzero  entries  because  b\, ... ,  share  the  same  relaxation  frequencies,  as  a 
consequence  of  the  orientation  and  location  invariance.  This  means  the  matrix  X  should 
be  row  sparse,  having  only  nonzero  entries  on  certain  rows. 

In  other  words,  the  invariance  property  of  the  relaxation  frequencies  translates  to  the 
row-sparse  property  of  the  matrix  X.  Row  sparsity  is  precisely  the  property  one  can  exploit 
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in  an  algorithm  that  estimates  the  DSRF  of  a  target  from  multiple  measurements. 

3.1.1  Row- Sparsity  Measure 


To  promote  solutions  that  are  row-sparse,  a  measure  is  needed  to  quantify  the  degree  of 
row  sparsity.  One  straightforward  measure  is  to  collapse  the  matrix  X  row-wise  by  the 
vector  norm  ||  •  ||g  [47],  resulting  in  a  sparse  vector.  This  can  be  done  via  the  function 
7 Z(q  :  Mm+1xL  — y  Mm+1  defined  by 


j(m)(X) 


Em=l  loS  (I  lxm| \q  +  e)  »  P  =  0 

ll^9WII?  =  ESl|Xm||?  ,0<p<l  ’ 


(33) 


where  e  >  0  is  a  small  positive  real  number  introduced  for  stability. 

The  purpose  of  the  log-sum  function  and  the  tp  quasi-norms  in  (33)  is  to  approximate 
the  t0  quasi-norm  (||  •  ||o),  which  is  the  number  of  nonzero  entries  in  a  vector.  Fig.  28  gives 
the  contour  plots  of  the  different  tp  quasi-norms  as  well  the  log-sum  function.  The  log-sum 
function  is  a  closer  approximation  to  the  to  quasi-norm  than  the  t\  norm  and  other  p  <  1 
quasi- norms.  The  log-sum  function  has  a  slope  that  vanishes  near  the  axes  like  ||  •  ||o  does. 
On  the  other  hand,  the  t\  norm  is  less  similar  to  ||  •  ||o-  Other  functions  have  also  been 
suggested  to  better  approximate  the  to  quasi-norm  [39,40,50]. 
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(c) 


(d) 


Figure  28:  Surface  plots  of  various  (quasi-)norms  in  M2.  (a)  £o  quasi-norm,  (b)  i\  norm, 
(c)  Elog(|^|  +  e).  (d)  fd/2  quasi-norm. 


While  there  are  conditions  that  guarantee  minimizing  the  i\  norm  is  equivalent  to  finding 
the  minimum  £o  quasi- norm  [51,52],  these  conditions  are  not  practical  for  DSRF  estimation 
when  there  is  measurement  noise  and/or  an  imperfect  model  (dictionary).  On  the  other 
hand,  it  can  be  shown  that  minimizing  the  £p  quasi-norm  for  p  <  1  or  the  log-sum  function 
via  reweighted  algorithms,  gives  sparser  solutions  and  more  robust  performance  [39,40]. 

3.1.2  MMV  Optimization 

Using  the  row-sparsity  measure  J^p’q\X)1  the  invariance  of  the  relaxation  frequencies  can 
be  exploited  to  effectively  model  an  EMI  signal.  To  estimate  the  model  parameters,  one 
optimizes  the  problem 

a.rgmin  \\A'X  -  B'\\f  +  A Jip’q){X),  (34) 


9te(A) 

91c  (B) 

where  A'  = 

,  B'  = 

3m(A) 

3m  (B) 
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1 1  •  1 1  p  is  the  Frobenius  norm  and  A  is  a  regularization  parameter  that  controls  the  trade-off 
between  sparsity  and  the  modeling  residual. 

This  optimization  (34)  for  MMV  is  similar  to  that  of  the  SMV  (Section  2.1.3)  but  with 
the  vectors  replaced  by  matrices,  the  vector  norm  replaced  by  the  Frobenius  norm,  and  the 
£p  norm  replaced  by  j(p’q\X).  The  real  and  imaginary  parts  are  again  separated  to  ensure 
real- valued  solutions. 

Equation  (34)  is  in  the  form  of  a  regularized  MMV  problem,  and  can  be  solved  by 
iteratively  reweighted  MMV  techniques,  such  as  M-FOCUSS  or  M-IRLl  [53,54].  Algorithms 
to  solve  this  problem  are  examined  in  Section  3.2,  and  their  performance  is  compared  in 
Section  3.3,  where  M-FOCUSS  is  chosen  as  the  default  solver  for  its  robustness. 

The  above  optimization  requires  selecting  the  regularization  parameter  A,  as  explained 
in  Section  2.5,  which  is  nonintuitive.  An  alternative  optimization  is  to  use  a  threshold, 
0  G  Rl,  to  specify  the  difference  between  the  modeled  signals  and  the  measurements,  which 
usually  contains  unwanted  signals  such  as  noise  or  soil  response.  The  alternative  threshold 
MMV  method  is 


argmin  J^p'q>(X)  subject  to  Co  (B  —  AX)  <6,  (35) 

x 

where  Cg  :  CNxL  — >  is  a  function  like  1Zgq  but  operates  along  the  columns  of  a  matrix. 

That  is, 

Ceq(E)  =  [||ei||?,...,||eL||?]T,  (36) 

/  N  \  l/q 

where  I \el\\q  =  f  \en,l\9j  •  (37) 

This  form  provides  more  physical  intuition  such  as  limiting  the  measurement  noise.  How¬ 
ever,  the  solution  to  (35)  is  more  sensitive  to  6  than  the  regularized  approach.  This  threshold 
MMV  can  be  iteratively  solved  similar  to  the  implementation  of  the  regularized  M-IRL1  [54]. 


3.1.3  Interpolation  of  Relaxation  Frequencies 


Sampling  in  the  relaxation  frequency  domain  always  introduces  the  issue  of  an  actual  re¬ 
laxation  Q  occurring  between  two  consecutive  sampled  relaxation  frequencies  Ca  and  Cb- 

Ca<  C  <  Cb-  (38) 
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Once  the  relaxation  frequencies  have  been  identified,  the  spectral  amplitudes  can  be  found 
via  a  standard  linear  least-squares  minimization  using  the  estimated  relaxation  frequencies 
(.  The  estimated  (interpolated)  relaxation  frequencies  are  denoted  by  £1,  (2,  ■  ■  ■ ,  where 
K  is  the  estimated  number  of  relaxation  frequencies.  Because  X  is  row-sparse,  K  <C  M . 
When  the  estimate  is  accurate,  K  ~  K,  if  not  K  =  K. 

First,  a  coefficient  matrix  Z  6  CNxK  is  constructed: 

^  1-.M/C1  1-.M/C2  i-jm/Cjc 

l+jwi/6  I+JW1/C2  '  '  '  1+j'wi/Ck- 

^  Wgg/Ci  i-j^  2/C2  2/Ck 

I+JW2/6  I+JW2/C2  '  '  '  1+^2/Ck 


^  l-j^jy/Cl  1— j^jy/C2  1—  JUn/Ch 

1+jwjv/Cl  I+JWJV/C2  l+JWjv/Cx 
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Unlike  the  dictionary  A  that  is  underdetermined,  Z  is  likely  overdetermined,  i.e.,  N  >  K. 
Therefore,  the  residual  from  the  linear  least-squares  minimization 

min  |  |h/ —  Zq|||  (42) 

Q 

is  unlikely  to  be  zero.  This  residual,  which  is  perpendicular  to  the  column  space  of  Z ,  is 
regarded  as  the  noise/unwanted  signal  rejected  by  the  EMI  model.  In  addition,  because 
the  estimated  (  are  sparse  and  adjacent  (  are  interpolated,  Cfc  are  spaced-out,  implying 
the  columns  of  Z  are  linearly  independent.  Because  Z  is  overdetermined  and  has  linearly 
independent  columns,  there  is  a  unique  solution  to  the  linear  least  squares  (42),  which  can 
be  obtained  via  the  pseudoinverse. 

To  ensure  only  real-valued  spectral  amplitudes,  a  real-valued  matrix  Z  £  K2Arx^  is 
constructed  similar  to  that  of  A'  in  (18).  The  same  arguments  for  Z  being  overdetermined 
and  having  linearly  independent  columns  transfer  to  Z  .  Therefore,  the  spectral  amplitudes 
can  be  estimated  using  the  pseudoinverse  Z  + : 

ct  =  Z+b’  l  =  1 . . .  L,  or  (43) 

C  =  Z  B'.  (44) 

For  each  measurement  b;,  the  estimated  DSRF  is  Si  =  {(C k,Ck,i)  ■  k  =  The 

estimated  relaxation  frequencies  are  /-independent,  invariant  from  measurement  to  mea¬ 
surement. 

3.2  MMV  Optimizers 

Several  algorithms  have  been  proposed  to  solve  the  MMV  problem.  In  this  section,  these 
solvers  are  examined  to  select  the  ones  that  are  best  suited  for  the  DSRF  estimation.  Five 
algorithms  are  considered:  M-BP  [47,55],  M-OMP  [47],  M-FOCUSS  [53],  ReMBo  (with 
BP)  [56],  and  M-IRL1  [54], 

A  set  of  noiseless  simulations  were  performed  to  evaluate  these  algorithms.  The  simu¬ 
lation  results  suggest  that  the  non-convex  optimizations,  M-FOCUSS  (p  =  0)  and  M-IRL1, 
deliver  the  most  robust  performance.  These  two  methods  are  further  examined.  Regulariza¬ 
tion  is  introduced  to  accommodate  noisy  measurements.  When  regularized,  the  simulation 
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result  using  noisy  measurements  suggest  that  M-FOCUSS  is  slightly  better  than  M-IRL1. 
In  addition,  because  M-FOCUSS  takes  less  computation  time,  it  is  chosen  to  be  the  default 
solver  for  DSRF  estimation  in  this  work. 

3.2.1  Algorithms 

Most  of  the  MMV  algorithms  are  extensions  of  existing  SMV  recovery  methods.  Extension 
of  the  Basis  Pursuit  to  the  MMV  problem  (M-BP)  is  considered  in  [47,55],  where  the 
objective  is  to  minimize  the  number  of  rows  containing  nonzero  entries  while  satisfying 
B  =  AX.  The  problem  is  formulated  as 

min  ||77^(X)||o  subject  to  B  =  AX ,  (45) 

where  ||  •  ||o  is  number  of  nonzero  entries  in  a  given  vector.  As  in  the  SMV  problem,  (45) 
is  NP-hard  but  can  be  made  convex  as  an  i\  minimization  problem: 

min  ||77^(A)||i  subject  to  B  =  AX .  (46) 

When  the  nonzero  rows  in  X  are  sparse  enough,  (46)  recovers  the  same  solution  as  (45). 
The  condition  for  which  (45)  and  (46)  are  equivalent  can  be  found  in  [47,55].  It  was  also 

shown  that  an  exact  recovery  does  not  depend  on  the  iq  norm  chosen  for  TZ^q. 

On  the  other  hand,  greedy  algorithms  have  also  been  extended  to  accommodate  MMV 
problems  [53,57,58].  Various  MMV  methods  based  on  Matching  Pursuit  (MP)  have  been 
proposed,  such  as  the  MMV  orthogonal  matching  pursuit  (M-OMP).  The  condition  for 
exact  recovery  was  also  established  [47,53]. 

From  a  slightly  different  approach,  the  ReMBo  method  proposed  in  [56]  solves  a  MMV 
problem  by  recasting  it  into  a  series  of  SMV  problems.  The  method  can  incorporate  both 
convex  relaxation  and  greedy  algorithms,  and  is  shown  to  be  robust. 

Sparsity  could  be  further  enhanced  through  iteratively  reweighting.  In  particular,  it 
was  shown  in  [39, 49]  that  sparse  solutions  for  a  SMV  problem  can  be  found  via  iterative 
reweighted  least-squares  (IRLS),  with  which  the  FOCUSS  algorithm  [59]  is  identified.  A 
M-FOCUSS  algorithm  that  extends  FOCUSS  to  the  MMV  problem  was  introduced  in  [49]. 
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The  M-FOCUSS  algorithm  solves  the  problem 


minJ^p,2)(X)  subject  to  B  =  AX.  (47) 

In  addition,  it  was  also  shown  in  [40]  that  sparse  solutions  can  be  obtained  via  an  iterative 
reweighted  t\  algorithm  (IRL1),  when  minimizing  a  log-sum  objective  function.  As  a  part 
of  this  research,  the  IRL1  algorithm  is  extended  to  the  MMV  case  (M-IRL1)  [54],  which 
iteratively  reweights  the  M-BP  algorithm.  The  M-IRL1  algorithm  solves  the  problem 

min  J^'l\x)  subject  to  B  =  AX .  (48) 

The  following  section  examines  the  performance  of  these  MMV  algorithms. 

3.2.2  Algorithm  Performance  —  Noiseless 

To  assess  the  performance  of  each  MMV  algorithm,  a  simulation  using  noiseless  measure¬ 
ments  is  conducted.  Noisy  simulations  are  considered  in  later  sections  that  follow.  Let 
N  =  20,  M  =  30,  and  L  =  5.  The  entries  of  the  real-valued  matrix  A  are  independent  and 
identically  distributed  standard  Gaussian  random  variables.  The  multiple- measurement 
vectors  are  constructed  by  B  =  AX o  where  Vo  has  K  rows  with  nonzero  entries.  The 
locations  of  the  K  rows  are  selected  uniformly  at  random,  and  the  nonzero  entries  of  Vo 
are  drawn  as  in  A.  An  exact  recovery  is  obtained  when  X  =  Xo-  The  above  simulation  is 
repeated  500  times  per  MMV  method.  From  the  simulation  result  (Fig.  29),  it  is  observed 
that  M-FOCUSS  ( p  =  0)  and  M-IRL1  have  the  highest  empirical  rate  of  exact  recovery  at 
different  K.  These  two  methods  are  further  examined  for  estimating  the  DSRF. 


Figure  29:  Performance  of  MMV  methods  for  noiseless  data. 
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3.2.3  M-IRL1  and  M-FOCUSS 

It  is  shown  here  for  noiseless  measurements,  the  special  case  of  M-FOCUSS  with  p  =  0 
and  M-IRLl  solve  almost  the  same  problem.  The  slight  difference  is  in  the  choice  of  the 
row-norm  7 Zg  ,  where  M-FOCUSS  uses  an  ( 2  row-norm  and  M-IRL1  uses  an  t\  row-norm. 
The  M-IRL1  algorithm  [54]  minimizes  the  log-sum  objective 

J((U)(X)  =  ^  log  (||xm||i  +  e)  (49) 

m=  1 

through  an  iteratively  reweighted  t\  algorithm.  The  formulation  of  the  algorithm  is  devel¬ 
oped  in  Appendix  B.  In  short,  the  algorithm  performs  a  series  of  t\  minimizations: 

X.^  =  argmin||W^7^£1(X)||i  subject  to  B  =  AX, 

x 

where  W1-^  =  diag[icf\  w^\  ■  ■  ■ ,  w^+1].  (50) 

The  weights  are  updated  by 

w™+1)  =  TT - (On  T  ’  m  =  1, . . . ,  M  +  1.  (51) 

||xmW||i  +  e 

On  the  other  hand,  the  M-FOCUSS  algorithm  [53]  minimizes  the  diversity  measure 

J^2\X)  =  \\nea(X)\\P  =  l|xm||2,  0  <  p  <  1  (52) 

m=  1 

subject  to  B  =  AX.  By  solving  the  Lagrangian  derived  from  (52),  an  iterative  reweighted 
least  squares  is  obtained  [53]: 

X(i)  =  W(<)AH(AW(i)AH)-1.B,  (53) 

where  =diag[w^, . . . ,  wj^j,  , ] ,  the  superscript  H  denotes  the  Hermitian  transpose,  and 
the  weights  are  updated  by 

^m+1)  =  (l|xmW||2)P  •  (54) 

When  p  =  0,  the  weight  is  (||xmW||2)  2,  and  M-FOCUSS  is,  in  effect,  minimizing 

j(°.2)(X)  =  ^log||xm||2,  (55) 

m=  1 
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which  is  very  similar  to  the  objective  of  M-IRL1,  J^0,1^(X).  This  is  true  because,  it  can  be 
readily  shown,  the  partial  derivative  of  J^0,2)(X)  with  respect  to  xmj  is 

/  r.  \  V2 


9J(°> 2)(X)  _  d 


dx 


m,l 


dx 


m’Z  rn= 1 


—  I  |xm|  1 2  %m,h 


J=  1 


and  when  substituted  into  the  derivation  of  M-FOCUSS  in  [53],  the  coefficient  of  xmj., 

1 1 xm 1 1 2^ 2 ,  is  identified  with  the  update  weight,  which  is  the  same  as  the  weight  in  (54)  with 
p  =  0.  That  is,  by  iterating  (53)  using  the  update  equation  (54)  with  p  =  0,  the  M-FOCUSS 
algorithm  effectively  minimizes  a  log-sum  objective. 

The  M-FOCUSS  algorithm  with  p  =  0  and  M-IRL1  are  almost  equivalent  except  that 
M-FOCUSS  uses  lZe2  (q  =  2),  whereas  M-IRLl  is  derived  using  1 However,  the  choice 
of  q  for  the  row-norm  should  not  make  a  significant  difference  because  the  log-sum  behaves 
similar  to  the  quasi-norm,  which  counts  the  number  of  nonzero  elements  regardless  of 
the  magnitude.  Indeed,  Chen  and  Huo  [47]  showed  that  under  certain  conditions,  the  row- 
norm  can  be  any  arbitrary  vector  norm,  and  the  problem  min  ||7?^  (-X’)Ho  (45)  can  be  solved 
exactly  by  solving  the  problem  min  ||77^(X)||i  (46). 

3.2.4  Regularization 

To  account  for  noise  in  the  measurements  and  modeling  error,  regularized  versions  of  M- 
FOCUSS  and  M-IR.L1  are  used.  To  regularize  M-IRL1,  the  iterated  minimization  with  an 
equality  constraint  in  (77)  is  replaced  with  the  following  relaxed  problem: 

=  argmin||AX-R||F  +  A||Ww77^(X)||i.  (56) 

x 

The  regularization  parameter  A  balances  the  emphasis  between  the  modeling  error  and  the 
sparsity. 

The  regularized  M-FOCUSS  iteratively  computes 

xW  =  wW  Ah(AW«  Ah  +  AI)_1R,  (57) 

where  I  G  WNxN  is  the  identity  matrix  [53]. 
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Since  regularization  effectively  loosens  the  equality  constraints,  the  regularized  algo¬ 
rithms  can  accommodate  actual  relaxation  frequencies  that  are  not  in  the  dictionary.  Sim¬ 
ulations  are  performed  to  assess  the  performance  of  each  regularized  algorithm  in  the  pres¬ 
ence  of  noise.  Comparisons  between  the  two  algorithms  are  shown  in  Section  3.3.3  along 
with  other  numerical  simulations. 

3.3  MMV  Simulations  with  Synthetic  Data 

Simulations  are  conducted  to  verify  the  validity  of  the  proposed  MMV  method  and  to 
gain  understanding  of  its  behavior  under  various  conditions.  Simulations  demonstrated 
here  include  a  four-relaxation  target  (Section  3.3.1),  the  performance  of  MMV  vs.  SMV 
(Section  3.3.2),  performance  comparison  between  different  solvers  (Section  3.3.3),  and  per¬ 
formance  comparison  for  various  number  of  measurements  L  (Section  3.3.4). 

In  Section  3.3.3,  M-FOCUSS  (p  =  0)  is  found  to  be  robust  and  fast  to  solve  the  regu¬ 
larized  MMV  optimization  problem  (34).  Thus,  for  the  rest  of  the  work,  unless  otherwise 
specified,  M-FOCUSS  (p  =  0)  is  assumed  to  be  the  MMV  solver.  The  regularization  pa¬ 
rameter  A  chosen  by  a  simple  function  of  SNR  is  described  in  Section  3.6. 

3.3.1  Synthetic  Four-relaxation  DSRF 

The  proposed  MMV  method  is  tested  using  a  four-relaxation  DSRF  synthesized  at  60  dB 
SNR  with  additive  white  Gaussian  noise.  The  spectral  amplitudes  are  assigned  according 
to  an  uniform  distribution  between  —1  and  1.  Ten  measurements  ( L  =  10)  are  used.  The 
estimation  result  using  M-FOCUSS  is  shown  in  Fig.  30.  The  relaxation  frequencies  £  as  well 
as  the  spectral  amplitudes  are  correctly  recovered.  All  samples  share  the  same  estimated 
C  because  of  the  sparsity-promoting  term  The  average  estimation  error  (EMD)  is 

0.014  decades. 

The  same  set  of  data  is  also  processed  using  one  measurement  per  estimate  (SMV).  The 
result  is  shown  in  Fig.  31  where  the  SMV  estimates  are  seen  to  have  more  variation  in  the 
estimated  £.  The  averaged  EMD  is  0.035  decades,  which  is  still  quite  low. 
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Figure  30:  Theoretical  and  MMV-estimated  DSRF  of  a  four-relaxation  target.  The  solid 
lines  indicate  the  true  relaxation  frequencies.  The  stems  with  solid  dots  are  the  actual 
spectra;  stems  with  hollow  circle  markers  are  the  estimates. 
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Figure  31:  Theoretical  and  estimated  DSRF  of  a  four-relaxation  target.  Blue  circles  are 
the  MMV  estimates;  red  squares  the  SMV  estimates. 

3.3.2  Performance  vs.  Signal  to  Noise  Ratio 

To  study  how  the  proposed  method  performs  in  noise,  a  Monte  Carlo  simulation  is  per¬ 
formed.  The  true  spectrum  is  from  a  target  with  a  four-relaxation  DSRF  including  nega¬ 
tive  relaxation  coefficients  Ck ■  A  set  of  L  measurements  is  synthesized  using  the  same  four 
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relaxation  frequencies  randomly  distributed  in  [wi,  w/v].  The  spectral  amplitudes  are  ran¬ 
domly  generated  but  constrained  to  have  the  same  sign  per  relaxation  frequency  (k-  This 
sign  constraint  improves  the  results  for  the  SMV  method,  because  the  L  measurements  are 
averaged  together  to  get  the  SMV  data  which  increases  its  SNR. 

Figure  32  shows  the  performance  of  the  MMV  and  SMV  methods  using  the  synthesized 
measurements.  The  simulation  is  set  up  with  L=10  measurements  per  trial  and  1000  trials 
per  SNR.  The  simulation  result  suggests  that  the  MMV  method  has  a  significant  advantage 
over  the  SMV  approach  even  though  the  simulation  is  setup  to  favor  the  SMV  method. 
Compared  to  the  performance  of  the  SMV  method,  the  MMV  gives  a  smaller  EMD  error 
by  about  half  a  decade  at  a  given  SNR. 


Figure  32:  Monte  Carlo  simulation  on  goodness  of  estimation  vs.  SNR.  The  error  bar 
indicates  the  10th  and  90th  percentile. 

3.3.3  Solver  Performance  Comparison 

A  simulation  is  conducted  to  examine  the  performance  of  two  algorithms  for  solving  the 
MMV  problem:  regularized  M-FOCUSS  and  M-IRL1.  The  simulation  is  similar  to  the  one 
in  Section  3.3.2  but  with  L  =  20  and  K  =  4.  The  Monte  Carlo  sample  size  is  1000  per 
SNR.  The  simulation  result  shows  that  M-FOCUSS  on  average  delivers  a  lower  estimation 
error  than  M-IRL1  (Fig.  33).  In  addition,  from  the  recorded  computation  time  (Fig.  34), 
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M-FOCUSS  runs  about  1000  times  faster  than  M-IRL1.  The  computation  time  of  both 
methods  is  observed  to  be  independent  of  the  SNR. 

The  M-IRL1  algorithm  is  implemented  using  the  CVX  package  which  implements  convex 
optimization  under  MATLAB  [60].  Because  CVX  implements  a  general  convex  optimization 
method,  further  speedup  for  M-IRL1  can  be  achieved  via  a  customized  implementation. 


Figure  33:  Monte  Carlo  simulation  on  goodness  of  estimation  vs.  SNR.  L  =  20  and  K  =  4. 
Sample  size  is  1000  per  SNR. 


Figure  34:  Monte  Carlo  simulation  on  computation  time  per  estimation  vs.  SNR.  L  =  20 
and  K  =  4.  Sample  size  is  1000  per  SNR. 


In  an  earlier  investigation  of  this  research  [54],  a  simulation  similar  to  the  one  in  Sec¬ 
tion  3.3.2  is  performed,  but  the  optimization  problem  is  complex-valued,  and  the  simulation 
uses  relaxation  frequencies  that  must  be  in  the  dictionary.  That  is,  for  any  actual  Ck  there 
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exists  a  Cm  such  that  (k  =  Cm-  Because  the  relaxation  frequencies  are  in  the  dictionary, 
it  is  possible  to  exactly  recover  them  even  under  noise.  In  this  case,  the  performance  is 
measured  by  the  exact  recovery  rate.  In  [54]  M-IRLl  is  observed  to  have  a  higher  probabil¬ 
ity  than  M-FOCUSS  to  exactly  recover  the  relaxation  frequencies.  This  is  because  in  [54], 
M-FOCUSS  is  not  constrained  to  consider  only  real-valued  answers,  which  degrades  its  per¬ 
formance.  On  the  other  hand,  M-IRL1  returns  only  real-valued  solutions  as  a  consequence 
of  its  implementation.  In  the  simulation  presented  here,  both  methods  are  constraint  to 
real- valued  solutions  by  separating  the  real  and  imaginary  part,  as  in  (34). 

3.3.4  Performance  vs.  Number  of  Measurements 

When  the  estimation  accuracy  versus  the  number  of  measurements  L  is  examined,  it  is 
observed  that  the  estimation  accuracy  increases  as  L  increases,  which  agrees  with  intuition 
(Fig.  35).  The  simulation  performed  here  is  similar  to  the  one  in  Section  3.3.2,  but  the 
optimal  regularization  parameter  is  used  instead  of  the  near-optimum  A  described  later  in 
Section  3.6.  The  intent  is  to  examine  the  best  possible  performance  for  a  given  L.  Because 
the  differences  in  performance  are  small,  the  optimal  A  is  used  to  help  differentiate  the 
performance  curves. 


Figure  35:  Monte  Carlo  simulation  on  estimation  accuracy  vs.  SNR  for  different  L,  as 
annotated.  The  number  of  relaxations  is  K  =  4.  Sample  size  is  1000  per  SNR. 
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3.3.5  Performance  vs.  p 


A  simulation  is  conducted  to  examine  the  estimation  accuracy  versus  SNR  for  different 
choices  of  p,  a  parameter  used  in  the  cost  function  j(p,q\X).  From  the  simulation  result 
(Fig.  36),  it  is  observed  that  there  is  a  significant  performance  difference  between  p  =  1 
and  p  <  1.  This  is  different  from  the  single- measurement  case,  where  the  performance  of 
p  <  1  is  somewhat  similar  to  that  of  p  =  1  (Fig.  27).  On  the  other  hand,  the  performance 
of  MMV  improves  significantly  once  p  is  less  than  1 . 

The  computation  time  versus  SNR  for  different  values  of  p  are  also  recorded  (Fig  37). 
It  is  interesting  to  note  that  the  computation  time  is  lower  for  lower  values  of  p. 


Figure  36:  Monte  Carlo  simulation  on  estimation  accuracy  vs.  SNR  for  different  p  values, 
as  annotated.  K  =  4.  Sample  size  is  1000  per  SNR. 


Figure  37:  Monte  Carlo  simulation  on  computation  time  per  estimation  vs.  SNR  for 
different  p  values. 
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3-4  Laboratory  Data 


The  functionality  of  the  MMV  estimation  method  is  verified  using  laboratory  data  where 
the  theoretical  DSRF  is  available.  An  automated,  non-metallic  measurement  facility  is  used 
to  measure  the  EMI  responses  of  a  target  at  various  positions  and  orientations  relative  to 
the  sensor  [44] .  The  proposed  MMV  method  is  observed  to  produce  accurate  estimates  and 
is  more  robust  than  the  single- measurement  counterpart. 

A  target  that  consists  of  three  mutually  orthogonal  copper  loops  is  examined.  The  loop 
diameters  and  thicknesses  are  3/20,  4/30,  and  5/36,  respectively  in  cm/ AWG (American 
Wire  Gauge).  The  target  is  configured  at  a  fixed  orientation  and  is  displaced  at  different 
positions  along  a  horizontal  axis  (x).  The  vertical  distance  between  the  target  and  the 
sensor  is  5  cm.  The  EMI  sensor  is  located  at  x  =  0.  The  response  of  the  target  is  measured 
at  nine  different  x  locations. 

For  the  MMV  method,  these  nine  responses  (L  =  9)  are  used  simultaneously  to  obtain 
an  estimate  of  the  DSRF.  For  the  SMV  method,  a  DSRF  is  estimated  at  each  location. 
Improving  the  SNR  by  averaging  the  nine  responses  is  not  as  effective  for  the  laboratory 
data  as  it  is  for  the  synthetic  data,  because  the  spectral  amplitudes  are  not  constrained  to 
be  of  the  same  sign  and  can  cancel  when  averaged,  lowering  the  SNR. 

The  estimates  from  these  two  methods  are  shown  in  Fig.  38.  Both  methods  produce 
satisfactory  estimates.  However,  the  MMV  method  gives  a  more  accurate  estimate  of  the 
relaxation  frequencies  than  the  SMV.  The  estimates  from  the  SMV  have  more  variation  in 
the  relaxation  frequencies  and  sometimes  extra  relaxations  are  introduced.  The  estimation 
error  of  the  two  methods  is  recorded  in  Table  1.  The  MMV  method  has  a  smaller  averaged 
EMD  of  0.127  decades  than  that  of  the  SMV,  0.157  decades. 


Table  1:  Estimation  error  (EMD)  of  the  three- loop  target 


x  (cm) 

-2 

-1.5 

-1 

-0.5 

0 

0.5 

1 

1.5 

2 

SMV 

0.1329 

0.0741 

0.3902 

0.3762 

0.1718 

0.0948 

0.0564 

0.0724 

0.0437 

MMV 

0.0750 

0.0461 

0.1716 

0.3231 

0.1915 

0.1231 

0.0886 

0.0704 

0.0505 
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Figure  38:  Estimation  of  the  DSRF  of  the  three  mutually  orthogonal  copper  loops  at  nine 
different  x  locations,  (a)  Black  diamonds  are  the  theoretical  DSRF,  blue  circles  the  MMV 
estimates,  and  red  squares  the  SMV  estimates,  (b)  A  3D  view  of  the  estimated  DSRFs  that 
shows  the  consistent  and  accurate  relaxation  frequencies  estimated  by  the  MMV. 
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3.5  Field  Data 


The  MMV  method  has  also  been  applied  to  the  field  data  where  measurements  are  corrupted 
by  various  factors,  including  the  soil  and  thermal  perturbations.  The  field  measurements 
used  in  this  section  are  the  same  set  of  data  used  in  Section  2.4,  where  targets  are  buried 
at  various  depths  in  grid  cells.  For  DSRF  estimation  using  multiple  measurements,  the  ten 
strongest  measurements  in  a  grid  cell  are  used  (L  =  10).  Consistent  DSRFs  are  observed 
for  targets  of  the  same  type,  indicating  the  functionality  of  the  proposed  MMV  method. 
Ten  types  of  targets  are  examined,  and  the  DSRFs  are  estimated  using  the  M-FOCUSS 
algorithm.  For  two  types  of  targets,  the  M-IR.L1  is  also  used  for  comparison. 

Type-A  and  Type-B  mines  are  modeled  using  both  M-FOCUSS  and  M-IRL1.  The 
description  of  the  mines  can  be  found  in  Section  2.4.  These  mines  have  moderate  to  strong 
EMI  responses.  The  results  of  the  estimated  DSRF  using  these  two  solvers  are  shown 
in  Fig.  39.  Both  methods  provide  stable  DSRF  estimates,  but  the  estimates  are  slightly 
different.  Because  the  actual  DSRF  is  unknown  for  the  field  targets,  it  is  not  clear  which 
solver  returns  more  accurate  estimates.  Nevertheless,  both  solvers  return  self-consistent 
DSRFs  that  can  be  used  as  a  feature  for  target  discrimination.  For  Type-A  mines,  the 
average  distance  among  the  estimates  are  0.047  decades  for  M-FOCUSS  and  0.079  decades 
for  M-IRL1.  For  Type-B  mines,  the  average  EMDs  are  0.088  decades  for  M-FOCUSS  and 
0.091  decades  for  M-IRLl. 

More  results  of  the  MMV  DSRF  estimation  are  shown  in  Fig.  40  and  Fig.  41.  As  with 
the  SMV  case,  the  DSRFs  obtained  using  the  MMV  method  are,  in  general,  consistent 
within  a  type.  Some  of  the  MMV-estimated  DSRFs  are  the  same  as  those  obtained  using 
the  SMV  method,  such  as  the  Type-W  mines  in  Fig.  18(d)  and  Fig.  40(d).  For  some  mines, 
such  as  the  Type-V  mines  in  Fig.  40(c),  the  MMV  method  seems  to  have  a  higher  resolution 
and  recovers  more  relaxations,  which  are  observed  in  the  SMV  only  when  the  SNR  is  higher. 
On  the  other  hand,  for  targets  like  the  Type-H  mines,  the  MMV  returns  DSRFs  that  are 
slightly  less  consistent  when  compared  to  the  SMV,  even  though  the  average  EMDs  are 
similar. 

From  the  field  data,  the  MMV  method  does  not  seem  to  show  a  significant  improvement 
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in  the  consistency  of  the  estimated  DSRF  when  compared  to  the  SMV  method.  This 
result  is  different  from  what  is  observed  in  the  synthetic  and  laboratory  data,  where  the 
MMV  achieves  a  higher  accuracy  in  the  estimated  DSRF  compared  to  the  SMV.  Overall, 
in  the  held  measurements  the  MMV  does  not  produce  more  consistent  DSRFs,  nor  does 
it  produce  less  consistent  DSRFs.  The  lack  of  the  expected  performance  gain  in  the  field 
measurements  might  be  caused  by  the  nonuniform  SNRs  across  the  measurements,  which 
affects  how  the  regularization  parameter  A  is  chosen.  In  the  following  section,  a  more 
sophisticated  A  selection  scheme  is  proposed.  For  this  section,  it  is  sufficient  to  demonstrate 
the  functionality  and  stability  of  the  MMV  DSRF  estimation  method. 
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Sample  Sample 


(c) 


(d) 


Figure  41:  Examples  of  estimated  DSRFs  from  field  targets  using  MMV.  The  average  EMD  listed  below  are  in  decades  and  SNR  in  dB. 
(a)  Eight  Type-C  mines;  EMD=0.19,  SNR=26.  (b)  Eight  Type-D  mines;  EMD=0.13,  SNR=82.  (c)  Eight  Type-E  mines;  EMD=0.12, 
SNR=26.  (d)  Four  Type-L  mines;  EMD=0.048,  SNR=85. 


3.6  Choosing  the  Regularization  Parameter  A 


A  simulation  is  performed  to  optimally  select  A  via  the  trade-off  between  the  modeling 
error  and  sparsity.  The  A-SNR  simulation  for  the  MMV  is  similar  to  that  of  the  SMV 
(Section  2.5),  except  that  there  is  an  additional  parameter  L  to  consider.  For  each  L,  the 
simulation  is  identical  to  the  SMV. 

As  in  the  SMV  case,  it  is  observed  that  the  optimal  A  that  gives  the  minimum  estimation 
error  is  quasi-independent  of  the  model  order  (Fig.  42).  In  addition,  the  EMD  (estimation 
error)  surfaces  are  well-behaved  and  smooth  with  respect  to  the  SNR  and  A.  To  obtain 
a  A-SNR  relation  independent  of  the  model  order,  an  average  of  the  EMD  surface  across 
different  model  orders  is  taken  (Fig.  43).  The  averaged  surface  is  smooth  and  has  a  “wide 
valley”  -  indicating  the  estimation  error  is  not  sensitive  near  the  optimal  A. 


Figure  42:  Monte  Carlo  simulation  of  the  goodness  of  estimation  (EMD)  of  spectra  of 
different  model  orders.  L  =  10. 

For  the  averaged  surface,  the  A  at  each  SNR  that  gives  the  minimum  EMD  is  roughly 
linear.  For  L  =  10,  a  linear  approximation  of  the  optimal  A  curve  is 

log  A  =  -0.059  •  SNR  -  2.44  (58) 

Better  approximations  can  be  achieved  through  higher  order  polynomials  or  splines. 
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Figure  43:  The  average  of  EMD  surfaces  of  various  model  orders,  i.e. ,  the  average  of  the 
surfaces  in  Fig.  42.  The  curve  with  asterisk  markers  traces  out  the  optimal  A’s.  The  line 
with  square  markers  approximates  the  optimal  A’s.  The  number  of  measurements  is  L  =  10. 

The  MMV  A  selection  is  in  fact  more  complicated  than  that  of  the  SMV  because  mea¬ 
surements  can  have  different  SNRs.  The  simulation  suggested  above,  however,  assumes  a 
uniform  SNR  across  measurements.  To  bridge  the  gap  between  the  uniform  SNR  simulation 
and  the  nonuniform  SNR  which  occurs  in  practice,  a  simple  mean-regularization-parameter 
A  can  be  used.  That  is, 


(59) 


i=i 


i=i 


where  A i  is  the  regularization  parameter  for  /th  measurement.  This  mean-A  works  well  when 


the  SNRs  do  not  vary  drastically,  and  it  is  demonstrated  in  the  field  data  the  functionality 
of  this  method. 

A  more  comprehensive  A-SNR  simulation  can  be  performed  to  account  for  the  varying 
SNRs,  which  may  lead  to  more  consistent  field  DSRF  estimates.  Two  factors  can  be  included 
in  the  current  simulation  setup.  The  first  is  to  assign  the  multiple  measurements  with  a 
nonuniform  distribution  of  signal  power  while  keeping  the  noise  power  constant.  A  reason¬ 
able  distribution  could  be  parabolic-shaped,  which  resembles  the  signal  power  distribution 
in  the  recorded  field  measurements  (e.g.,  Fig.  15).  The  second  factor  is  the  soil  response  in 
the  simulation.  Either  synthesized  or  field-measured  soil  responses  can  be  considered. 
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CHAPTER  IV 


TARGET  DISCRIMINATION 


The  concept  of  discriminating  targets  based  on  the  relaxations  of  a  target  is  well  known  [5, 
12],  but  the  inability  to  obtain  good  DSRF  estimates  has  been  a  limiting  factor  [30].  Using 
the  DSRF  estimation  methods  developed  in  this  work,  it  is  now  possible  to  discriminate 
targets  based  on  the  DSRF. 

Recall  that  the  EMI  response  of  a  target  can  be  modeled  by  the  expression 

ff  (<*)  =  «,  +  !>  ,  (60) 

where  the  relaxation  frequencies  C  are  intrinsic  to  the  target  and  are  invariant  with  respect 
to  orientation  and  location.  The  relaxation  frequencies  naturally  serve  as  a  good  signature 
or  fingerprint  for  a  target.  The  spectral  amplitudes  Ck  change  with  respect  to  orientation 
and  location,  so  the  are  usually  not  a  ideal  signature.  However,  in  the  case  of  landmine 
targets,  most  landmines  are  buried  at  fixed  orientations  so  the  changes  in  are  consistent. 

In  this  work,  the  strategy  is  to  use  both  the  relaxation  frequencies  and  the  spectral 
amplitudes  for  target  discrimination,  i.e. ,  use  the  DSRF  S  =  {(OoQ,)  :  k  =  1 . . .  K}.  The 
methods  developed  in  this  chapter  can  be  readily  adjusted  to  discriminate  targets  using 
only  the  relaxation  frequencies  by  replacing  the  spectral  amplitudes  with  a  constant.  This 
constant  should  be  normalized  such  that  Y!k=ick  =  1- 

To  visualize  how  target  discrimination  based  on  the  DSRF  is  possible,  consider  the  plots 
shown  in  Fig.  44  for  several  instances  from  two  types  of  landmine.  It  is  clear  the  two  mine 
types  have  very  different  DSRF  fingerprints  -  Mine  Type-A  has  two  relaxations  and  Mine 
Type-B  has  six  relaxations.  In  addition,  the  spectral  amplitudes  for  each  type  of  mine  are 
self-consistent.  Therefore,  it  is  possible  to  differentiate  the  types  of  targets  based  on  the 
DSRF.  More  examples  of  targets  exhibiting  different  DSRFs  can  be  found  in  Section  2.4  and 
3.5.  While  multi-class  classification  is  possible,  in  this  chapter  only  two  classes,  landmine 
and  non-landmine,  are  considered  to  give  a  working  example  of  DSRF-based  discrimination. 
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Figure  44:  Estimated  DSRF  of  real  landmines.  The  spectral  amplitude  is  represented  by 
the  intensity:  darker  the  gray,  larger  the  amplitude,  (a)  Eight  Type-A  mines:  low  metal 
content,  nonmagnetic,  and  moderate  EMI  response  antipersonnel  mines,  (b)  Eight  Type-B 
mines:  medium  metal  content,  magnetic,  and  strong  EMI  response  antipersonnel  mines. 


To  quantify  the  dissimilarity  between  two  spectra,  the  Earth  Mover’s  Distance  (EMD) 
is  used.  For  example,  the  EMD  can  be  used  to  quantify  the  distance  between  Type-A  and 
Type-B  mines  by  averaging  the  distance  between  all  pairs  of  DSRF  from  the  two  types. 
This  can  be  done  for  instances  from  all  types  of  mine  available  as  well  as  clutter  objects. 
By  doing  so,  a  similarity  map  is  obtained  (Fig.  45).  The  similarity  map  is  a  symmetric 
matrix  of  EMD  between  all  pairs  of  DSRFs.  The  diagonal  is  zero  because  that  is  the  EMD 
between  a  DSRF  and  itself  which  is  zero. 

The  map  shows  that,  in  general,  mines  of  the  same  type  are  similar  to  each  other, 
as  indicated  by  the  dark  blocks  on  the  diagonal.  On  the  other  hand,  targets  of  different 
types  are  dissimilar,  indicated  by  the  white  color  off-the-diagonal.  More  importantly,  land¬ 
mines  are,  overall,  different  from  the  clutter,  and  this  suggests  that  landmine  and  clutter 
discrimination  is  possible  using  the  DSRF. 
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Figure  45:  The  EMD  between  samples  from  about  one  hundred  independently  measured 
field  landmines  of  eleven  types  and  various  metal  clutter.  Darker  colors  denote  smaller 
distances  which  indicate  that  two  samples  are  more  similar. 


To  discriminate  targets  using  the  DSRF,  the  k-nearest  neighbor  (kNN)  algorithm  and 
the  support  vector  machine  (SVM)  are  considered  to  classify  targets.  The  EMD  plays  a  key 
role  in  formulating  these  DSRF  classifiers.  These  classifiers  are  presented  in  Section  4.1.1. 

To  aid  real-time,  on-the-field,  target  discrimination,  a  model-based  soil  prescreener  is 
introduced  to  identify  the  presence  of  metallic  objects.  The  DSRF  estimation  and  classi¬ 
fication  are  only  performed  when  a  target  is  present.  The  soil  prescreener  is  discussed  in 
Section  2.4.2. 

The  classifier  and  soil  prescreener  are  incorporated  in  a  landmine  detection  framework 
developed  in  Section  4.1.3.  Lastly,  the  performance  of  the  framework,  utilizing  the  DSRF 
classifiers  and  the  prescreener,  is  examined  and  compared  to  existing  detection  techniques 
in  Section  4.2. 
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4-1  Detection  Methods 

4.1.1  Target  Classifiers 

Two  DSRF-based  target  discrimination  techniques  using  the  kNN  and  the  SVM  are  devel¬ 
oped  here.  With  only  two  classes,  landmine  and  not-landmine,  the  classification  problem 
reduces  to  a  detection  problem. 

The  kNN  predicts  the  class  of  a  target  based  on  the  labels  of  its  closest  k  neighbors. 
In  the  case  of  the  kNN,  the  distance  measure  used  to  quantify  the  similarity  between  two 
DSRFs  is  the  EMD.  The  Euclidean  distance  is  a  poor  measure  of  the  similarity  between 
the  two  DSRF  because  the  relaxations  of  two  DSRFs  are  usually  sparse  and  not  aligned,  as 
discussed  in  Section  2.2.1. 

In  the  case  of  the  SVM,  given  a  target’s  DSRF  parameter  set  S,  the  target  is  classi¬ 
fied/labeled  using  the  decision  function 

f(S)  =  sign  ^  aiytK(S 'T,  S)  +  ^  ,  (61) 

where  K  is  the  kernel  (explained  shortly),  Stp  the  training  data,  y*  G  {  —  1,  +1}  the  training 
class  labels,  at  the  trained  weights,  and  f3  the  trained  threshold.  Only  a  few  a’s  are  nonzero, 
i.e. ,  the  a’s  are  sparse.  The  Sip  that  correspond  to  the  nonzero  a*  are  called  the  support 
vectors. 

For  the  kernel,  a  generalized  radial  basis  function  based  on  the  EMD  is  used  [61]: 

KEmd(Si,S2)  =  exp(— pEMD(Si,  S2)),  (62) 

where  p  is  a  scaling  parameter.  For  brevity,  (62)  is  called  the  EMD  kernel  [62],  While 
it  has  not  been  possible  to  prove  that  the  EMD  kernel  satisfies  Mercer’s  condition,  i.e., 
Kemd  is  positive  semi-definite,  it  is  observed  that  the  EMD  kernel  is  positive  semi-definite 
in  practice  [62].  Furthermore,  it  should  be  noted  that  kernels  that  do  not  satisfy  Mercer’s 
condition  can  still  perform  well  [61]. 

4.1.2  Soil  Prescreener 

A  soil  prescreener  based  on  the  soil  model  (21)  is  presented  here.  The  prescreener  filters  out 
responses  that  are  like  those  due  to  the  magnetic  properties  of  the  soil.  This  reduces  the 
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computation  cost  in  unnecessary  DSRF  estimation.  The  prescreening  process  is  performed 
very  efficiently  via  a  linear  least-squares  minimization. 

In  Section  2.4.2,  a  soil  model  is  proposed  based  on  the  log-linear  trend  in  the  soil 
frequency  response: 

Hg(u)  =pi  +P2  (hicu  +  j^)  ,  (63) 

where  p\  and  P2  are  model  parameters.  The  real  part  has  a  linear  trend  with  respect  to 
the  log-frequency,  and  the  imaginary  part  is  a  constant.  Recall  that  for  a  soil  response  bo 
measured  at  N  frequencies, 

bo  =  Gp  +  noise  ,  (64) 


where 

1  Incur  +  j  7t/2 
1  lncu2  +  j  7t/2 

and  p  = 

1  lncujv+j7r/2 

Using  the  system  of  equations  (64),  a  measured  soil  response  can  be  fitted  efficiently  via  a 
least-squares  minimization: 

bG, m  =  G(GuG)-1GnbG.  (65) 

The  magnitude  and  residual  of  soil  responses  for  6000  samples  collected  at  locations 
that  are  reported  to  have  no  metal  content  are  shown  in  Fig.  46.  Some  magnitudes  are 
strong  (>  — 125  dB)  because  metal  targets  were  actually  present  nearby.  For  most  samples, 
the  magnitude  of  the  modeling  residual  is  noticeably  smaller  than  the  response,  indicating 
that  the  samples  fit  the  soil  model  well.  Ground  responses  that  do  not  fit  well,  i.e.,  have 
high  residuals,  often  occur  when  strong  metallic  objects  are  near  by. 

Because  the  model  describes  a  behavior  very  specific  to  the  soil,  the  model  can  be  used 
as  a  prescreener  to  determine  whether  a  target  is  present  based  on  the  fitting  residual. 
Given  a  frequency  response,  the  response  is  fitted  to  (63)  using  (65)  and  the  fitting  residual 
e  is  compared  with  a  threshold  6  to  determine  whether  a  metallic  object  is  present: 

e  true  if  e  >  9 

target  present  =  <  (66) 

false  otherwise. 
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A  reasonable  choice  of  9  for  our  measurement  is  — 135  dB,  as  suggested  by  Fig.  46. 
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Figure  46:  Samples  of  blank  responses  fitted  to  the  soil  model  (21).  The  samples  are 
ordered  so  that  response  decreases  with  increasing  sample  number. 

4.1.3  Detection  Framework 

A  detection  framework  that  incorporates  the  soil  prescreener  and  the  target  classifier  is 
presented  here.  The  framework  is  designed  to  be  suitable  for  practical  applications  where 
measurements  are  acquired  sequentially  in  real-time.  It  is  designed  according  to  the  scenario 
where  a  detection  vehicle  carrying  the  EMI  sensors  is  driven  forward  and  EMI  responses  bn 
are  collected  sequentially. 

The  prescreener  first  screens  out  the  responses  that  are  absent  of  metallic  objects.  Re¬ 
sponses  that  pass  the  prescreener,  indicating  a  target  is  present,  are  then  processed  to 
estimate  their  DSRFs.  Based  on  the  estimated  DSRF,  the  classifier  then  labels  the  re¬ 
sponses  as  landmine  or  not-landmine.  The  use  of  the  prescreener  significantly  reduces  the 
amount  of  data  processed  by  the  DSRF  estimator  and  the  target  classifier.  Because  the 
prescreener  takes  very  little  computation  time  compared  to  the  estimator  and  the  classifier, 
the  average  computation  time  is  also  greatly  reduced  by  using  the  prescreener. 

A  simple  voting  mechanism  is  employed  to  discourage  temporary  mislabeling  of  land¬ 
mines  by  taking  advantage  of  the  sequential  measurements.  As  the  detection  vehicle  passes 


i 

■ 

< 


0  1000  2000  3000  4000  5000  6000 

sample 


77 


over  the  target,  often  multiple  measurement  are  collected  consecutively  for  that  target.  In 
this  case,  multiple  labels  Xi  are  produced,  and  a  more  confident  decision  can  be  made  based 
on  the  recent  labels.  A  landmine  is  determined  to  be  present  only  when  p  out  of  the  past  q 
labels  are  marked  as  landmine.  The  voting  rule  reduces  the  false-alarm  rate  and  increases 
the  confidence  level. 

The  proposed  framework  is  summarized  below  and  a  flow  chart  is  shown  in  Fig  47. 


Detection  Framework _ 

Input:  bn,  6  ,  p,  xn-q+i . . .  xn-\ 

Output:  decisionn 

1  Fit  bn  to  soil  model  (21)  and  obtain  residual  e. 

2  if  e  <  6  then 

3  |  xn  =  0 

4  else 

5  S  =  estimated  DSRF  of  bn 

6  xn  =  classify(S)  (0  or  1) 

7  if  E”=n-g+ 1  xi  >  P  then 

8  |  decision^  =  1 

9  else 

10  decisionn  =  0 

11  return  decisionra 


Figure  47:  Flow  chart  of  the  detection  framework. 
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4  ■  2  Performance 


The  proposed  method  is  applied  to  a  data  set  acquired  from  a  testing  field  that  contains 
26  types  of  landmines  and  various  types  of  metallic  and  nonmetallic  clutter.  The  field  is 
divided  into  220  grid  cells  where  targets  are  buried.  This  is  the  same  testing  field  described  in 
Section  2.4.  About  145  EMI  responses  are  collected  per  grid  cell.  In  total,  32,148  responses 
are  collected  for  the  whole  field.  The  acquisition  hardware  used  is  described  in  [8]. 

The  EMI  responses  are  collected  sequentially  as  the  detection  vehicle  is  driven  down 
the  lane,  and  the  responses  are  fed  into  the  detection  framework  described  earlier.  The 
parameters  are  chosen  such  that  p  =  10,  q  =  20,  and  0  =  — 135  dB.  A  snapshot  of  the 
output  of  the  framework  is  shown  in  Fig.  48.  Because  the  responses  are  filtered  [8],  a  target 
response  has  multiple  lobes, e.g.,  grid  71  to  74. 
o 

-100 
-200 

Figure  48:  A  snapshot  of  the  output  of  the  detection  framework.  The  labels  on  the 
horizontal  axis  indicate  the  grid  numbers,  and  the  vertical  axis  is  the  response  magnitude 
in  dB.  The  curved  lines  are  the  strength  of  the  responses  measured  at  21  frequencies.  The 
target  types  are  noted  above  the  grid  number.  The  blue  lines  (near  -200  dB)  indicate  points 
that  are  marked  as  soil;  black  dots  (at  -100  dB)  indicate  points  that  are  labeled  as  landmines; 
red  diamonds  indicate  a  declaration  of  landmine. 

In  Fig.  48,  it  is  seen  that  the  soil  prescreener  is  quite  effective  and  the  voting  rule  reduces 
false  alarms.  In  grid  72,  a  miscellaneous  clutter  is  labeled  as  a  landmine  for  a  few  times, 
but  because  the  number  of  mislabels  is  small  (<  p),  a  landmine  is  not  declared  and  a  false 
alarm  is  avoided. 

More  output  of  the  detection  framework  is  shown  in  Fig.  49,  where  multiple  lanes  are 
displayed.  In  general,  landmines  are  correctly  detected,  labeled  by  the  red  diamonds,  the 
soil-only  responses  are  correctly  identified,  indicated  by  the  blue  lines,  and  the  metallic 
clutters  are  correctly  classified.  In  grid  66,  the  target  is  weak  and  is  screened  out  by  the 
soil  prescreener.  Upon  close  examination,  this  target  has  a  soil-like  frequency  response. 
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Figure  49:  Output  of  the  detection  framework  when  applied  to  multiple  lanes.  The  notation  is  the  same  as  in  Fig.  48. 


The  performance  of  the  detection  framework  using  the  SVM  is  summarized  in  Fig.  50. 
The  receiver  operating  characteristic  (ROC)  curve  achieves  a  high  detection  rate  of  0.96  at  a 
low  false-alarm  rate  of  0.10.  Other  operating  points  also  provide  satisfactory  performances. 


Figure  50:  ROC  curves  of  the  proposed  method  and  that  of  Fails  et  al..  For  the  kNN 
ROC  curve,  k  =  7. 

The  detection  process  does  not  take  much  computer  time.  With  a  single  pre-trained 
SVM,  the  whole  test  field  (32,148  responses)  can  be  classified  using  the  above  process, 
including  estimating  the  DSRF,  in  30  seconds  on  a  2.66  GHz  CPU  with  960  MB  RAM. 

The  performance  of  the  detection  framework  using  the  kNN  is  also  shown  in  Fig.  50. 
The  performance  is  comparable  to  that  of  the  SVM.  However,  the  processing  time  is  much 
longer  (10+  minutes)  due  to  the  many  distance  computations  required  to  find  the  nearest 
neighbor  per  measurement.  While  the  kNN  may  not  be  suitable  for  real-time  application, 
it  is  quite  robust  when  sufficient  training  data  are  available,  as  is  the  case  here.  When 
training  data  are  scarce,  the  SVM  is  likely  to  have  smaller  generalization  error. 

The  classifiers  are  trained  per  grid  cell  using  a  leave-one-out  cross-validation  (LOOCV), 
i.e. ,  the  classifiers  are  trained  at  each  grid  cell  using  responses  from  other  grids  cells.  Only 
the  strongest  responses  in  a  grid  cell  are  used  for  training.  Grid  cells  containing  only  soil 
are  excluded  in  the  training  process. 

The  performances  of  the  DSRF-SVM  and  DSRF-kNN  are  compared  to  the  detection 
performance  of  a  simple  energy  detector,  also  shown  in  Fig.  50.  The  energy  detector  does 
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not  require  training  and  declares  a  landmine  present  if  the  energy  of  the  response  is  greater 
than  a  threshold.  To  avoid  bias  in  the  real  part  of  the  response,  which  is  due  to  the  DC 
magnetization  of  the  soil,  the  energy  is  computed  using  the  mean  of  the  imaginary  part 
of  the  response.  Compared  to  other  methods,  the  energy-detection  performance  is  rather 
poor. 

Figure  50  also  includes  the  ROC  curve  for  the  method  of  Fails  et  al.  [6]  where  the 
performance  is  evaluated  on  the  same  data  set  using  LOOCV.  While  the  proposed  method  is 
slightly  better,  the  simulation  done  by  Fails  et  al.  does  not  utilize  sequential  measurements. 

All  the  ROC  curves  in  Fig.  50  saturate  at  0.99  detection  rate.  This  is  due  to  the 
misclassification  of  one  particular  landmine.  Upon  close  examination,  it  is  found  that  the 
response  of  this  landmine  is  very  weak  and  is  indistinguishable  from  the  soil  response.  In 
the  proposed  framework,  this  target  is  filtered  out  by  the  soil  prescreener.  Other  features, 
such  as  the  soil  response  and  the  magnetic  property  of  targets,  may  be  included  to  provide 
even  more  robust  performance.  The  study  presented  here  only  demonstrates  the  strong 
potential  of  using  the  DSRF  for  landmine  detection. 
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CHAPTER  V 


CONCLUSIONS 


In  this  thesis,  two  effective  modeling  techniques  for  EMI  signals  are  developed:  one  for  a 
single  measurement  vector  and  one  for  multiple  measurement  vectors.  These  techniques 
overcome  several  long-standing  obstacles  in  estimating  the  discrete  spectrum  of  an  EMI 
response.  This  is  achieved  by  first  linearizing  the  problem  through  enumeration  and  then 
exploiting  the  notion  of  sparsity.  The  multiple-measurement  technique  exploits  the  invari¬ 
ance  of  the  relaxation  frequencies  to  further  improve  the  estimation  performance.  The  two 
techniques  provide  robust  features  for  a  DSRF-based  classification  that  was  developed  to 
demonstrate  the  strong  potential  of  DSRF-based  subsurface  target  detection. 

The  main  contribution  of  this  research  is  a  new  technique  in  estimating  the  DSRF  of 
a  target  from  its  EMI  response,  which  enables  DSRF-based  subsurface  target  detection. 
Discriminating  targets  based  on  DSRF-like  models  has  been  studied  for  decades,  but  the 
inability  to  estimate  the  DSRF  reliably  has  been  the  primary  limiting  factor  until  advances 
made  by  this  research.  The  developed  methods  fall  in  the  domain  of  detecting  and  discrim¬ 
inating  subsurface  targets  using  broadband  EMI  sensors.  The  contribution  of  this  research 
has  been  to  first  characterize  the  EMI  response  of  a  target  through  modeling  and  then 
discriminate  targets  based  on  the  characterizations  using  classification  algorithms. 

In  Chapter  1,  the  motivation  and  background  of  subsurface  target  discrimination  using 
broadband  EMI  sensors  are  provided.  Two  classes  of  broadband  EMI  model  are  introduced 
-  continuous  and  discrete.  It  is  argued  that  the  discrete  model  is  valuable  for  target  discrim¬ 
ination  due  to  the  position  and  orientation  invariance  of  the  discrete  relaxation  frequencies. 
In  this  chapter,  the  difficulties  in  extracting  the  parameters  of  the  discrete  model,  i.e. ,  es¬ 
timating  the  DSRF,  are  addressed.  Previous  works  and  their  shortcomings  are  examined. 
A  similarity  map  based  on  the  DSRF  for  different  types  of  targets  is  given  to  demonstrate 
the  potential  to  classify,  hence  discriminate,  targets  based  on  the  discrete  model. 
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In  Chapter  2,  a  sparsity-regularized  DSRF  estimation  method  using  a  single  measure¬ 
ment  is  developed.  The  estimation  problem  is  first  linearized  through  enumeration  in  the  re¬ 
laxation  frequency  domain.  The  linearized  problem  is  then  solved  by  a  sparsity-regularized 
least-squares.  Through  tests  performed  using  synthetic,  laboratory,  and  field  data,  this 
method  is  found  to  give  satisfactory  estimates.  A  soil  model  is  proposed  based  on  observed 
soil  responses,  and  an  augmented  dictionary  that  includes  the  DSRF  and  the  soil  model  is 
developed  to  model  the  target  and  soil  response  simultaneously.  A  simulation  is  performed 
to  study  the  regularization  parameter  A  and  then  derive  a  simple  formula  for  selecting  A  as 
function  of  SNR. 

In  Chapter  3,  the  DSRF  estimation  method  proposed  in  Chapter  2  is  generalized  to  ac¬ 
commodate  multiple  measurements.  Exploiting  the  invariance  of  the  relaxation  frequencies, 
a  matrix  equation  is  formulated  with  the  property  of  row-sparsity.  Using  this  row-sparsity,  a 
robust  DSRF  estimation  method  is  developed  for  the  multiple-measurement  case.  Through 
simulations  and  tests  using  the  laboratory  and  field  data,  this  proposed  method  is  found  to 
be  robust  and  fast.  A  performance  gain  over  the  single  measurement  method  is  observed 
over  a  wide  range  of  SNRs  in  the  simulated  results. 

In  Chapter  4,  the  DSRF  is  applied  to  discriminate  targets  using  classification  algorithms 
in  conjunction  with  the  EMD.  A  soil  prescreener  is  introduced  to  filter  out  responses  that 
lack  metallic  objects.  A  detection  framework  is  developed  to  include  the  soil  prescreener, 
the  DSRF  estimation  method,  and  the  classifier.  A  sliding-window  voting  mechanism  is  in¬ 
troduced  to  discourage  temporary  mislabeling.  The  DSRF-based  target  detection  is  demon¬ 
strated  to  deliver  high  detection  but  low  false-alarm  rate.  The  detection  framework  provides 
a  working  example  for  DSRF-based  subsurface  target  detection. 

Two  ideas  for  future  research  are  suggested  in  Section  2.4.2  and  Section  3.6.  The  first  is 
to  investigate  the  effect  of  soil  and  to  quantify  the  range  of  signal-to-soil  ratio  in  which  the 
DSRF  estimation  performs  satisfactory.  The  second  is  to  perform  a  more  comprehensive 
regularization  parameter  simulation  that  accounts  for  the  nonuniform  SNRs  in  multiple 
measurements.  The  results  of  this  future  research  can  further  improve  the  robustness  of  the 
DSRF  estimation  methods. 
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APPENDIX  A 


THE  EMI  MODEL 


The  physical  electromagnetic  model  of  the  EMI  response  and  its  relation  to  the  discrete 
EMI  model  is  presented  here.  It  is  also  discussed  here  the  relationship  among  several  forms 
of  the  EMI  model,  including  time-domain  and  frequency-domain  model. 

A.l  A  Physical  Model  for  the  EMI  Response 

The  EMI  response  of  a  metallic  object  is  the  result  of  the  interaction  between  the  trans¬ 
mitting  loops,  the  receiving  loops,  and  the  object’s  magnetic  polarizability  M  [63].  It  can 
be  shown  by  using  reciprocity  that  response  H (u>)  is 

H(u)  =  aHRTM(w)HT,  (67) 


where  a  is  a  real  constant,  Ht  is  the  magnetic  held  generated  by  the  transmitting  loop,  Hr 
is  the  magnetic  held  of  the  receiving  loop  if  it  is  driven,  and  M  is  a  complex,  frequency 
independent,  second-rank  tensor. 

Equation  (67)  can  be  expanded  because  the  magnetic  polarizability  of  a  target  can  be 
written  as  a  sum  of  relaxations  [64]: 


M(uj)  =  T0T0-^Tk( 


k= 1 


ju/Ck  ' 
1  +  jA/Cfc' 


(68) 


where  T).  is  a  real  constant  and  T*.  is  a  real,  symmetric,  second-rank  tensor.  The  hrst  term 
is  due  to  the  bulk  magnetic  permeability  of  the  target,  which  is  assumed  to  be  frequency 
independent,  and  the  second  term  is  due  to  the  currents  induced  in  the  target  [63]. 
Expanding  M  in  (67)  using  (68),  the  response  becomes 

JA/Cfc  v 


H(u)  =  «Hrj 


T0To-^!7Mr 

k=i 


+  ju/Ck' 


Ht 


=  a 


T1  U  T rj- 1  u  \  '  rp  (  3^  / Ck  Tf|i  u 

JotlR  1  QjlT  ~  /  .  !/-  1 


k= 1 


1  +  Jw/ Ck ' 
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(69) 


do  —  ^  dk 


k= 1 


ju/Ck 
1  +  ju/Ck 


where  dk  =  oTk Hr7  T^Ht-  Note  that  is  independent  of  position  and  orientation  while  dk 
is  dependent  of  the  position  and  orientation  of  the  target  relative  to  the  sensors.  This  form 
provides  more  intuition  about  the  physical  process  of  EMI,  where,  again,  the  first  term  is 
due  to  the  bulk  magnetic  permeability  of  the  target  and  the  second  term  is  due  to  the  eddy 
currents  induced  in  the  target. 

From  (69),  a  unit-step  time  response  can  be  obtained  using  the  inverse  Laplace  trans¬ 
form: 

h(t)  =  C~ 

=  C~ 

=  dou(t)  —  dke~tc,ku[t)  (70) 

k=X 


A. 2  Forms  of  the  EMI  Model 


The  EMI  model  (69)  can  be  written  in  other  forms  that  are  useful  for  signal  processing. 
Two  alternative  forms  are  examined  in  this  section.  When  the  nonnegative  least-squares 
method  was  first  developed  in  this  research,  the  following  form  was  introduced  [33]: 


where  the  parameters  are 

The  relationship  between 
H(u)  =  do 
=  do 


H(u)  =r0  +  Y l 


Tk 


related  to  the  physical  form  (69)  by 

do  =  r0  +  y^yrk  and  dk  =  rk. 
k=  1 

(69)  and  (71)  is  as  follows: 
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(72) 
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In  the  sparsity-promoting  methods,  another  form  (6)  was  introduced  to  give  uniform- 
normed  columns  in  the  linearized  model,  which  is  related  to  (71)  by 


H(u)  =  ro  +  Y  y 


?'k 
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.  Cfc  rfc 

where  c0  =  r0  +  2^  y  and  ck  =  — . 


k= 1 


(73) 

(74) 


This  uniform-norm  model  (73)  is  relate  to  the  physical  form  (69)  using  (72)  and  (74): 


do  =  co  +  Y  ck  and  dk  =  2 ck.  (75) 

k= 1 

The  form  in  (73)  is  preferred  for  the  DSRF  estimation  because  has  a  unit 

norm  independent  of  Cfc,  so  the  dictionary  in  (16),  populated  by  enumerating  ,  has 

columns  with  equal  norms,  implying  that  each  Cm  is  equally  likely  to  be  selected  by  the 
DSRF  estimator.  This  is  true  because  the  proposed  DSRF  estimation  (18)  minimizes  the 
£p  norm,  and  the  tp  norm  discourages  large  entries.  If  a  column  has  a  smaller  norm,  then 
its  corresponding  entry  in  the  weighted  selection  vector  needs  to  take  on  a  larger  value  to 
compensate  the  small  column  norm,  which  is  penalized.  Therefore,  selecting  a  column  with 
a  smaller  norm  is  penalized  more  compared  to  other  columns  with  larger  norms.  For  this 
reason,  a  dictionary  with  uniform  column  norms  is  desired  and  the  unit-norm  form  for  H (uj) 
(73)  is  preferred. 
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APPENDIX  B 


JOINTLY  SPARSE  VECTOR  RECOVERY  VIA  ITERATIVELY 
REWEIGHTED  ix  MINIMIZATION 

The  sparse  minimization  problem  (45)  can  be  approximated  by  minimizing  the  log-sum 

objective  function  on  the  norm  of  the  rows: 

M+ 1 

argminY^  log (||xm||g  +  e)  s.t.  B  =  AX ,  (76) 

* 

where  e  >  0  is  a  small  positive  real  number  introduced  for  stability.  Recall  that  ||xm||g  >  0 
are  the  entries  of  1Ziq(X),  the  iq- norm  of  the  rows  of  X  defined  in  (32).  For  simplicity, 
q  =  1  is  considered  here. 

Following  the  reweighting  scheme  in  [40],  (76)  can  be  solved  with  an  iterative  algorithm 
referred  here  as  M-IRL1: 

Step  1)  Initialize  count  i  =  0  and  win  =  1,  i  =  1, . . . ,  M  +  1. 

Step  2)  Solve  the  weighted  l\  minimization  problem 

XW  =  argmin  ||W^^77^1(A)||i  subject  to  B  =  AX, 

x 

where  =  diag[icf\  w%\  . . . ,  w^+1\.  (77) 

Step  3)  Update  the  weights: 

w^+1)  =  - ^ - ,  m  =  1, . . . ,  M  +  1.  (78) 

1 1  xm  Hi  T  e 

Step  4)  Terminate  on  convergence  or  when  i  reaches  a  specified  maximum  number 
of  iterations  imax.  Otherwise,  iterate  from  Step  2. 

While  (76)  better  promotes  sparsity,  it  is  nonconvex  and  a  unique  solution  is  not  guar¬ 
anteed.  The  proposed  M-IRL1  method  can  be  trapped  in  local  minima.  However,  when 
an  initial  point  is  properly  chosen,  the  algorithm  does  converge  to  the  global  minimum,  as 
shown  empirically  in  Section  3.3.  The  proposed  method  converges  as  argued  in  [50]. 
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B.l  Algorithm  Justification 


This  section  provides  a  justification  for  the  proposed  method.  A  similar  derivation  is  found 
in  [40],  which  is  the  following  with  L  =  1. 

In  (76),  substitute  for  ||xm||g  using  (31): 


arg  min  E  log  E  |a;m,i|+e]  s.t.  B  =  AX . 


(79) 


m= 1  \l= 1 


The  minimization  in  (79)  is  equivalent  to 


arg  min  E  log  E  umj  +  6  )  s.t.  s  G  C, 


(80) 


m=  1  \  1=1 


where  s  =  ( X ,  U)  and  C  is  the  convex  set  {(xmti,  umfi  j  B  =  AX  and  \xmj\  <  umj\. 

Recognizing  the  objective  function  in  (80)  is  concave,  which  is  below  its  tangent,  a  guess 
s G  C  can  be  improved  by  minimizing  a  linearized  objective  function  around  s^: 


s(fe+i)  =argmin^(s(fc))  -i-  Xg(s^)(s  —  s^)  s.t.  s  €  C, 

S 

where  g(s)  =  ^  log  ^  umj  +  e 

m= 1  \  1=1 

It  can  be  readily  shown  that 


From  (81)  and  (82), 


(X(fe+1),  £/(fe+1))  =  argmin  £ 

m= 1  Yll= 1  um,l  £ 

s.t.  B  =  AX  and  \xmj\  <  umjh 


(81) 


(82) 


(83) 


which  is  equivalent  to 


=  argmin  y  ^=1  s.t.  B  =  AX 

m=  1  <LjI= 1  \xm,l\  '  t 


Using  (31), 


X(fc+1)  =  argmin 


x- 


ill  1 1  q 


\(k) 


s.t.  B  =  AX. 


+  e 


(84) 


The  iteration  and  weights  in  (84)  define  the  proposed  algorithm. 
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APPENDIX  C 


EARTH  MOVER’S  DISTANCE 


Given  two  distributions  S  =  {(C i,(k)  ■  i  =  1 . . .  K}  and  S  =  {(C j,cj )  :  3  =  1  ...K},  the 
Earth  Mover’s  Distance  (EMD)  between  the  two  distributions  can  be  computed  by  solving 
the  optimization  problem  [43]: 


EMD  (S',  S)  =  rnir: 

sr^K  sr^K 
2^i=i  2^7=1 

L - - - - - 

fij  dij 

(85) 

fij 

sr^K  sr^K 
i=l  2~^j= 

■  1  fij 

K 

subject  to 

y  ~  g 

j= i 

i  = 

1... 

K 

(86) 

k 

y!  fv  ~  °j 

i= 1 

3  = 

1... 

.  K 

(87) 

K  I< 

k 

I< 

= 

=  min 

(E 

cn  cj  ) 

(88) 

i=l  j= 1 

i= 1 

j= 1 

fij  >0  i  = 

=  1 . . . 

K,j 

=  1 ...  A 

(89) 

where  fij  is  an  intermediate  variable  used  during  the  optimization  and  dij  is  the  distance 
function.  Adapting  the  illustration  in  Section  2.2.1,  S  is  the  piles  of  earth  and  S  the  holes. 
Equation  (86)  guarantees  no  overdraw  from  each  pile  of  earth,  (87)  guarantees  no  over  fill  at 
each  hole,  (88)  sets  the  problem  to  fill  up  the  holes  with  as  much  earth  as  possible,  and  (89) 
allows  only  moving  earth  into  holes  and  not  the  reverse. 

In  this  work,  spectra  should  be  normalized  having  sum  of  all  spectral  amplitudes  be 
unity  Q =  1).  In  this  case,  the  above  optimization  problem  is  simplified  to  having 
the  denominator  in  (85)  be  one  and  the  right-hand-side  of  (88)  be  unity.  The  EMD  also 
becomes  symmetric. 
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